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DEDICATION 


A Ginette 


ABSTRACT 


A new spherically symmetric numerical model for hailstone 
icing has been devised whichsimulates certain aspects of the discrete 
nature of the accretion process, and which considers the heat transfers, 
including internal conduction, to be time-dependent. As a first approx- 
imation, a uniform layer of given thickness is accreted instantaneously 
over the entire surface of a spherical hailstone with initially homo- 
geneous properties... Ihe subsequent history of the deposit: before the 
accretion of the next layer is divided into two stages. During the 
freezing stage, the deposit warms to OG owing to the formation of 
ice dendrites and it exchanges heat with the environment and the in- 
terior of the stone while the remainder of the supercooled water 
freezes. If the entire deposit is frozen before the accretion of the 
next layer, it enters the cooling stage, where exchanges of heat with 
the environment and the interior give rise to a changing deposit 
temperature. During both of these stages, the interior temperature 
profile of the hailstone is also calculated as a function of time. 

By making the hypothesis that the individual droplets accreted on 

the hailstone would have a temperature cycle identical to that of a 
uniform deposit, a simple estimate of the temperature distribution 

over the hailstone surface is obtained. From such a distribution, 

the mean surface temperature, the standard deviation and the proportion 


of the stone surface in wet growth are caiculated. 
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INTRODUCTION 


Hail formation can be studied through the analysis of hail- 
stones. However, in order to explain the very complex structure (both 
inside and at the surface) of hailstones, hypotheses and assumptions 
have to be made about the microphysical processes taking place during 
their initiation and growth inside thunderclouds. These assumptions 
and hypotheses lead to the formulation of conceptual models. 

Several microphysical models have been proposed (Schumann, 
1938; Ludlam, 1950; List, 1963) to simulate heat transfers during 
hailstone growth. These were able to explain some of the observed 
features of hailstone structure such as opaque ice, clear ice or spongy 
ice. The purpose of this work is to review these studies and to present 
an improved numerical model which simulates more realistically some 


aspects of hailstone icing. 
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CHAPTER 1 


EXISTING MODELS 


1.1 An Early Model 


Schumann (1938) was the first one to recognize that the 
surface temperature of an ice deposit formed by the accretion of 
supercooled water droplets should be warmer than the environment. 

He was also the first to consider that the latent heat of fusion re- 
leased on Freezing of the supercooled water had to be disposed of in 
one way or another. An energy budget was then established in which 
the latent heat released was transferred to the ambient air through 
conduction and convection and through evaporation of some of the 
water on the hailstone surface (sublimation of ice is not considered). 
In doing so the amount of heat conducted into the hailstone was 
neglected. By setting up a balance equation for heat transfers, 
whicn depended on the hailstone surface temperature, the cloud water 
content, the air temperature and the pressure, he was able to solve 
for an equilibrium surface temperature. Interestingly, his result 
indicated no dependence of the equilibrium temperature on the hail- 
stone dimensions. So for a given cloud model, the hailstone tempera- 


ture would be a function of height alone. This peculiarity is basically 
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due to an erroneous formulation of the ventilation coefficients used 
in expressing the heat and water vapour transfers. 

The basic assumptions of Schumann's model are that the hail- 
stone is a perfectly smooth sphere, and that it collects all the 
cloud water (but no ice particles) encountered in its path. Using 
results published by Bilham and Relf (1937), the terminal velocity of 
the stone is calculated to obtain the rate of accretion of supercooled 
water and therefore the rate at which latent heat is released, assuming 
total freezing. This should be equal to the rate at which the heat is 
dissipated through conduction to the air and through evaporation of 
water. In order to calculate the dissipative terms, Schumann postulates 
the existence of a ''stagnant layer of air'' around the falling hailstone 
inside of which al! transports are due to molecular processes. Accord- 
ing to his calculations the rate of heat dissipation per unit area has 
the same dependence on the radius as the terminal velocity so thas 
this variable drops out of the balance equation. 

Another important concept which this model gave rise to is 
the concept of critical liquid water content. Schumann was the first 
one to realize that there was a limit to the amount of supercooled 
water which could be frozen on the hailstone surface. Above that 
amount the rate of dissipation would be insufficient to balance the 
rate of production of latent heat. The hailstone would then warm up 
to 0°C and the excess collected water would not freeze. This limit 
is referred to as the critical liquid water content. Schumann 
thought that the excess collected water would be thrown off in the 


form of droplets. Again the value of the critical liquid water content 
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shows no dependence on the hailstone dimensions and is a function 
only of the cloud parameters. 

Although this model leads to erroneous conclusions, it was 
a first attempt to simulate realistically the heat transfers and it 
had the merit of leading to the very important concept of critical 


liquid water content. 


1.2 Present Models 


Using more accurate measurements of the drag and ventilation 
coefficients obtained with the aid of wind tunnels, Ludlam (1950, 1951, 
1958) corrected the heat balance equation set up earlier by Schumann. 
Essentially the same equation was used except for improved formulations 
of the heat transfers due to convection and conduction and due to 
evaporation. This results in an equilibrium temperature which is a 
Function of the hailstone size, as well as of the external cloud con- 
ditions. The critical liquid water content above which all the 
collected supercooled water cannot be frozen on the hailstone surface 
is also derived from this model, and it has a dependence on the hail- 
stone radius. When a hailstone enters a region for which the liquid 
water content exceeds the critical value, it enters a growth regime 
which Ludlam describes as the wet growth regime. In these conditions 
a film of liquid water covers the hailstone surface. Ludlam believed 
that most of this water would be shed by the hailstone. During this 
stage clear ice forms. Opaque ice occurs during the dry growth regime, 
where the cloud water content is subcritical and the accreted water 


freezes completely. The alternating layers of clear and opaque ice 
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characteristic of large hailstones are then explained by ''... successive 
passages through critical conditions caused by changes in radius, 
falling speed, temperature, pressure and liquid water content in the 
eloud=s! {Ludiam, 1950). 

This type of continuous model was given in its most complete 
form by List et al. (1963, 1965, 1967). The basic assumption of the 
model which is also an implied assumption of the earlier models, is 
the ate steady-state assumption. It is supposed that the rate of 
change of the hailstone radius is sufficiently slow that the heat 
transfer rates vary only very slowly and the hailstone remains at the 
equilibrium temperature. Again the heat transferred to the hailstone 
itself is assumed to be negligible. Calculations are done for spheri- 
cal hailstones, although later extended to spheroidal hailstones 
(1967). 

Under these assumptions, the genera! heat balance equation 
is set up as follows: 

OFF + Qe ee be Sw) (el) 
where Qe is the rate of heat exchange by conduction and convection; 


is the rate of heat exchange due to evaporation or sublimation; 


25 


Qep is the rate of heating of the accreted cloud droplets (in order 
for them to rise to the surface temperature); and Q. is the rate of 


release of latent heat by the freezing of the supercooled water drop- 


lets. 
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The originality of this ue besides the thoroughness in 
the expression of each term of the balance equation, is that it allows 
for incomplete freezing of the deposit. It is now known that spongy 
ice (a mixture of ice and water) could play an important role in the 
formation of hailstones. This model, assuming that none or only 
part of the non-frozen water is shed by the hailstone, permits one 
to make some conclusion about the conditions of formation of this 
particular type of ice. 

Upon writing each term as a function of the surface tempera- 
ture and of the ice fraction of the deposit, the balance equation 
can be rewritten as: 


! bs val es 
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ee , = On 
k: thermal conductivity of the air (cal cm 1 *c 1) 
eas O 
deposit equilibrium temperature (C) 
; fe) 
t,: air temperature (~C) 

Cc : constant which contains the product of mass transfer 
coefficient for water vapour and latent heat of 
vaporization or sublimation. It takes the value 
.207 cai °c cm 3 mb ! for liquid-gas transitions 
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intel 2. 26hep eit) ae cm > mb ! for solid-gas transitions. 
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T_: absolute temperature of the air (K) 

Sous oy! saturation vapor pressure over the hailstone 
surface and for the environment 

v: kinematic viscosity (cm2 s !) 

D: diameter of spherical hailstone (cm) 

8: roughness factor 

Ee §icollecti onwefifiicl ancy 

pi cloud water content (g cm 3) 

Le: latent heat of fusion of water (cal g !) 

c : specific heat capacity averaged from te to ty 

(callng ts) 


I: ice fraction of deposit 


For given cloud conditions, equation (1.2) relates D, 


Tet respectively, the hailstone diameter, the ice fraction, and 


q: 
the temperature of the deposit. Ignoring the finite freezing time 
Grothe deposit, Le= loonly if ty = c°c. Therefore for a given 
hailstone diameter and given cloud conditions, the temperature and 
the ice fraction (assuming no water is shed) are uniquely determined. 
Calculations were done for a ''typical'' cloud model relating air 
temperature, pressure and height. The results are presented in 
graphical form in the well-known ''Listograms'' (see Appendix 2). 
This model was also extended by List et al. (1967) to 
include melting. In that case the warmer environment would be 
giving heat to the stone, through conduction and convection, while 


condensation of water vapor on the colder hailstone surface wouid 


release latent heat of vaporization; this heat would then be used 
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to melt some of the ice on the hailstone. 

As mentioned earlier, the basic assumption of the models 
presented above is the quasi steady-state assumption. However, when 
a hailstone grows, its terminal velocity increases, which in turn 
increases its rate of accretion. *Since the rate.of accretion, and, 
therefore, the rate of latent heat release, increases faster than the 
dissipating terms, the surface temperature increases with time. 
According to the steady state models, this change in temperature 
should be sufficiently slow to be of no real significance. Never- 
theless, the effects of a rising surface temperature were investi- 
gated by Pellet and Dennis (1974) using a modified version of the 
continuous model. In this model heat transfers with the underlying 
hailstone are allowed. Writing the balance equation for this model, 


we have 
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where Me, C. 


ie are, respectively, the mass, specific heat 
capacity, and absolute temperature of the hailstone. Calculations 
were done for a hailstone growing in a stationary one-dimensional 
cloud model. Comparing the results with those obtained when 

Q; = 0, as in the steady-state models, it is found that there is 
very little difference between the two and that the stone adjusts 
fairly quickly to the equilibrium temperature. However, this last 


model assumes that the hailstone has a uniform temperature, and 


that the heat is conducted instantaneously into the stone interior, 
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which implies an infinite thermal conductivity. 


1237 AP Discrete Model 


The underlying assumption of all the models described 
above is that the accretion of supercooled water is continuous 
(whence the label continuous model). They treat the cloud water 
as if it were a gas whose density is the liquid water content. 

In reality the accretion process is discontinuous with the water 
being accreted in lumps as the individual droplets impinge on the 
hailstone surface. 

Macklin and Payne (1967) made a theoretical investiga- 
tion of ice accretion emphasizing the freezing of the individual 
droplets. Once a supercooled water droplet impinges on the surface 
of a hailstone, the equilibrium is disturbed and ice is nucleated. 
How rapidly and in which manner the ice phase will develop inside 
the droplet is not clear. It should depend greatly on the 
molecular mechanisms involved in the solidification process; but 
it should also depend on the rate at which the latent heat of fusion 
can be dissipated to the surroundings. Several authors (Macklin 
and Ryan, 1965) have investigated the growth of ice in bulk 
supercooled water at strong supercooling; in this case, because 
of the heat limitations, the growth takes a dendritic form. In 
the case of lower supercoolings (less than 2.5°C), the growth 
rate is slower and more orderly. Estimates of growth velocities 
of ice dendrites in supercooled water as a function of total 


supercooling are given in Figure l.1. 
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Figure 1.1 Linear crystallization velocities (LCV) as a function 


of supercooling, after Hobbs (1974). LCV are obtained 
by measuring the growth of ice in supercooled water 
contained in narrow tubes. The continuous line rep~ 
resents data obtained by Tammann and Blichner and the 
dashed line, data from Pruppacher. 
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We are not aware of eeeurarente of the growth rate of ice 
inside supercooled cloud droplets, but on the basis of Figure 1.1 one 
can expect that the initial growth should be very rapid. In order 
to study the freezing process in the dry growth regime of a hailstone, 
Macklin and Payne formulated the following model. 


"A newly arrived droplet (on the hailstone surface) 
is nucleated and its temperature rises rapidly towards 
0°C.... This will be called the initial freezing. 

. The subsequent freezing rate is controlled partly 
by heat conduction into the deposit and partly by 
forced convection and evaporation to the environment. 
When the freezing process is complete the droplet 
begins to cool by forced convection and in a steady 
Stave Itecools untily on thesaverage, the initial 
mean temperature of the deposit is regained before 
the arrival of the next droplet when the procedure 
is repeated. We distinguish therefore three main 
phases, involving different physical processes, 
during which the heat of a droplet is removed: the 
initial freezing, the subsequent freezing and the 
cooling phase. Although there is some overlapping 
of these three phases, the total time, Tt, it takes 
to remove the heat liberated by the droplet may be 
conveniently regarded as the sum of three individual 
times, namely, 


In the steady state situation t is then the average 
interval of time between the arrival of successive 
droplets." 

(Macklin and Payne, 1967, p. 201-202). 

In order to simplify the calculations, the authors consider 
the freezing of a uniform layer of supercooled water of thickness 6R 
on a sphere of radius R, initiaily at a uniform temperature Ty sO 
that heat flows will be in the radial direction only. Right at the 


outset it is assumed that t. is negligible with respect to te (al- 


though it is found later that this is not quite true for very thin 
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layers, less than 1 um) and that no heat is exchanged outside the 
deposit during that time. So in practice, the process starts with 
an ice-water deposit initially at 0°, and during the subsequent 
freezing, heat is conducted inside the sphere, and heat is trans- 
ferred to the air by forced convection and evaporation. The deposit 
is said to be frozen when all the latent heat available has been 


expended. Therefore, 
Heating -iha=s a0 20R ‘c ace (Te - 1.) § (135) 


The L.H.S. of equation (1.5) is the sum of the heat given 
to the environment at the rate H during time te and qe is the heat 
conducted inside the hailstone. The R.H.S. of equation (1.5) is 
the latent heat available after warming the accreted deposit to the 
melting point Us Macklin and Payne assumed without much evidence 


that Ht, is small enough to be negligible with respect to des also, 


£ 
in the evaluation of p> they neglected the spherical effects by 

using the solution for linear heat flow into a semi-infinite body 

having a plane surface. We will see later that this last assumption 

is not fully justified. Nevertheless, with the aid of these assumptions, 


Macklin and Payne obtain the following expression for the (subsequent) 


freezing time: 


2 
D 2 & 
t = mp2er2 {i ee | T) 
Aa a NSE A ie a EE (1.6) 
ip.c, v2 K 


where p7 is the density of water, 6R the deposit thickness, Le the 
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latent heat of fusion for water, C the specific heat capacity for 
water (although Le and care temperature dependent, this is not 


considered by the authors), V is the absolute value of T expressed 


d 
in Celsius, and Pir Ce, K are, respectively, the density, the specific 
heat capacity and the thermal conductivity for ice; 2 and ie are the 
air temperature and the melting temperature. 

Next comes the cooling phase. The authors' ideas about 
this stage are somewhat confusing. It seems as if they assume that 
the final deposit (hailstone) temperature Ty is determined by external 
conditions (and can be obtained by equations similar to those of the 
List model) and that the duration of the cooling phase ty must be just 
what is needed for the deposit to return to Ty3 the arrival of the 
next layer must then await the completion of this phase and the 
achievement of Toe In other words, the authors are investigating 
the processes by which a mean equilibrium deposit temperature is 
achieved, and they do not question this deposit temperature itself. 
The hypothesis is: In order to return the sphere to its initial state 
(i.e., where the average temperature is Ty)> the heat conducted into 
the sphere during the freezing process together with that given out 


by the ice layer in cooling from ie. to Ty have to be removed. Thus 


the heat per unit area to be removed is 


q. = 4 + e,,5Re, (7, = 1). (ler) 


The rate at which heat is dissipated to the environment depends on 
the difference between the surface temperature of the sphere and that 


of the air, as well as on the difference between the saturated vapor 
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pressure over the hailstone and that in the environment. Even though 
the dependence of the vapor pressure on the temperature is far from 
linear, the authors assume linearity, in order to be able to use the 
solution of the heat equation with the radiation boundary condition. 
No justification is given to support this assumption; it seems very 
difficult to evaluate the type of error introduced by this procedure. 
(Further discussion of this procedure is given in section 3.1.) 
Interestingly, results show that in order to remove all the 
heat absorbed during the freezing phase, the surface temperature 
falls below the original temperature Ty: Furthermore, it is found 
that this heat is dissipated some time before the arrival of the 
next layer so that the hailstone becomes somewhat colder than ini- 
tially. Table 1.1 gives a comparison of the time necessary to freeze 
the deposit and for the hailstone to return to the initial heat con- 
tent, and the time between the arrival of two successive layers. It 
is seen that the difference is quite significant especially for thick 
layers at cold conditions. This could mean that large cioud droplets 
could form deposits that are colder than those resulting from an 
equivalent number of smaller droplets. Although these findings were 
acknowledged by Macklin and Payne, no discussions of the reasons or 


the possible implications are presented. 
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Table 1.1 


Tt time necessary to freeze the deposit and for the 
hailstone to return to the initial heat content, 
compared to the time between the arrival of two 
successive layers, for two different thicknesses 

and different initial temperatures. t is the sum 

of t¢ and t. as calculated by Macklin and Payne 
(1967, Table 6). At is the time necessary for the 
hailstone to grow by the given thicknesses, using 
Ludlam's (1958, p. 49) formula. Tq is the initial 
temperature of the hailstone (1.5 cm in diameter) 
and the values correspond to the equilibrium tempera- 
ture calculated for the -20°C level of the Wokingham 
storm (Browning and Ludlam, 1962) for liquid water 
contents also given by Macklin and Payne (1967, 
Table 6). 
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1.4 Other Work with Internal Conduction 


The effects of internal conduction were also investigated 
by Hitschfeld and Stauder (1967). They calculated the temperature 
change of a hailstone falling in clear air and found that there was 
a considerable lag accompanied by a significant internal temperature 
gradient (2-3°C from surface to center). For a hailstone falling 
through supercooled rain, they found that the effects of internal 
conduction were over-shadowed by the latent heat effect of super- 
cooled water. 

Goyer et al. (1969) used the heat diffusion equation to 
calculate the rate of freezing of spongy hail. Their theoretical 
calculations were compared with experimental measurements showing 
very good agreement between the two. Applying their numerical model 
to a simplified thunderstorm model, they concluded that it is unlikely 
that sufficient time is available for large, spongy hailstones to 


freeze completely. 
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CHAPTER 2 


A NEW MODEL 


2.1 The Proposed Model 


All the continuous models described in Chapter 1 neglected 
the effects of a finite thermal conductivity for the hailstone. 
Macklin's results indicated the possibility of significant differences 
when the discrete nature of the accretion process was taken into 
account; however, only one accretion cycle was examined and the approx- 
imations made were not totally satisfactory. For these reasons, it 
was decided to develop a model which would further investigate the 
discrete nature of the freezing process by studying the freezing and 
cooling cycles of several accreted layers with the proper heat trans- 
fers. This model would also look into the effects of allowing heat 
conduction within the hailstone due to internal temperature gradients. 

The proposed model is basically an extension of the model 
developed by Macklin and Payne, and described in the first chapter. 
The two major differences are that the accretion cycle is repeated 
successively several times, with the temperature profile resulting 
from the previous cycle being the initial condition for the next one. 


Also, the heat conduction equation, together with the boundary conditions 
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is solved by means of numerical fees, thereby eliminating the 
necessity of making approximations to the analytic solutions in order 
to make them readily solvable. 

The main assumptions of the model are spherical symmetry 
(spherical hailstone), and constant density and thermal conductivity 
throughout the hailstone interior. It is also assumed that the 
deposit has the same properties as the hailstone, and that it is 
accreted uniformly over the surface. Because of these assumptions 
the mathematical problem is spherically symmetric: the temperature 
and ice fraction are functions of the radius and time and we are 
concerned with radial heat transfers only. 

| The model considers a spherical hailstone growing under 
constant cloud conditions. Melting and refreezing are not investi- 
gated; however, the case of incomplete freezing in the wet growth 
regime is taken into account. 

Under the continuous growth hypothesis, the time necessary 
for a spherical hailstone to grow by the radial increment 6R is given 


by (Ludlam, 1958) 


oR , Ay) 


where V, is the terminal velocity of the stone, We the liquid water 
ConLcntemeathe collection craiciency. P; the density of the deposit. 
In our model it is assumed that a deposit of thickness 6R is in- 
stantaneously accreted; therefore, t is taken to be the mean time 
between the arrival of each layer, or the time duration of one cycle; 


alternatively this could be interpreted to be the mean time between 
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the arrival of two droplets over the same area on the hailstone 
surface. 

The choice of Pas the density of ice, for the deposit 
density is justified by the fact we are assuming the deposit to have 
the same properties as the hailstone. It is not at all clear what 
should be the density of a hailstone since it could vary greatly, 
depending on the surface temperature at the time of accretion 
(Macklin, 1962). At any rate, the value of 917 kg m 2 was selected, 
as suggested by List (1963), although values as low as 300 kg m 2 
have been measured for rime ice formed on cylinders under very cold 
conditions (Macklin, 1967). 

For a hailstone, the collection efficiency can be set 
equal to unity (Ludlam, 1958; List, 1963). However this ignores 
possible water shedding at high liquid water contents in the wet 
growth regime (List et al., 1976). Substituting in equation (2.1) 


we get 


oR 
T= 9g. 000 vp (Si0e (22) 


When a hailstone is at the terminal velocity, the drag 
force exerted by the air is equal to the weight of the hailstone. 
For a falling sphere the terminal velocity is given by (is tectea Loe 


1967) 


g p,D 2 
Veen bes iS (ST jae (2.3) 
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hailstone density, D the hailstone diameter, om the air density, Cy 
the drag coefficient. For hailstones of diameter larger than 5 x 10 3 
m, the drag coefficient is nearly constant and can be set equal to 
0.50 (List et al., 1967). However, Strong (1974) suggested a value 
closer to 0.6. 

Substituting in equation (2.2) we get 


4 
m2 are ere (a (2.4) 
p 


Here, T is the air temperature in degrees Kelvin, R is the hailstone 


radius in m, and p is the pressure in kPa. 


2.2 The Accretion Cycle 


Within one cycle, three processes are taking place. First, 
the deposit is accreted instantaneously. It is assumed that dendrites 
are nucleated quasi-instantaneously and uniformly within the deposit 
so that a negligible amount of heat is exchanged outside the deposit. 
The validity of this assumption has been discussed by Macklin (1967) 
and was analyzed in the first chapter. Therefore, all the latent heat 
released at that stage is used up in warming the deposit to 0°c. Once 
the deposit has reached the melting temperature, no more freezing is 
possible without some mechanism for dissipating the released latent 


heat and the resulting ice fraction is given by 


veto pe At sae ee (2.5) 
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where is is the air temperature, Se and a are respectively the water 
and the ice specific heat capacity averaged from te to 0°C. This 
formula supposes that the latent heat released on freezing of the 
fraction I of the supercooled water at the temperature te is used up 
in warming the entire deposit, ice and water alike, to o°c. It is 


equivalent to the following shorter expression 


ore ite 
a —_—s (2.6) 
L,(0 ¢)) 
However the physical meaning here is not realistic. It 


seems in fact to imply that all the supercooled water warms to O16. 
and that the heat needed for this warming is provided by freezing a 
certain portion I of the deposit at 0°C which is impossible. This 
initial instantaneous phase corresponds to the initial freezing stage 
of Macklin's model. 

Next comes the freezing phase, analogous to Macklin's 
subsequent freezing. During this period, the remaining liquid water 
(or part of it) will be frozen. The total mass of unfrozen water is 


given by 


M. = AnR?6Rp, (1 - I) , . (2.6) 


and the total available latent heat is given by 


Oo 
Ho = Lp(o°c)M. = Ank2srp, (1 - I) Le(O'C) . (2-7) 


20 


a ae amie . 
a ids ae 


7 Wy ia Ps at eee? = hd : 

- . 7 : 

7 i} 7 ® 

ene oo: 2 wire) .: 4 Pea hinweed )eek wef 
try a 


re 
>a q 1 a P <¥e 7 : vin > O35 mt oe 


= ‘ - ' 
~ = » 
: " ner 
- 
Li : 7 a 4 i — 
' - 5 : S . ~ 
| : 7 ; ?' P os owe 
<* : “ : . ; EBs a 
1] _ 
v7 ‘VTS nari ¢ Ppt “es ‘we Heh weg, 
“i. ae 
woth | wipe? tel heat 
{q@ lat Sd ope 
7 : cee a rs. : 
vet 4 . s , en a 
% > ' 7 ' a a 6 i : 7 7 
A ; iad : - ’ t F Gs 
bg 7 
- : : 7 a : : et ed 
le lhe, RA os 
7 i ’ a te: ea n a cs 


ea 7 
ties rm lig rel aia ita 7 


1, i ee 
Pale bes ng mr 
‘ aan a % Vs — 


2] 


Here again p; is used according to the initial assumption 
that the deposit has the same density as the hailstone. The deposit 
temperature is maintained at 0°C as long as there is liquid water 
present: the available latent heat is disposed of by exchanges 
with the colder environment and with the hailstone interior. If all 
the deposit becomes frozen before the arrival of the next layer, that 
is, if all the latent heat available has been released, the hailstone 
enters the cooling phase. The name itself is quite explicit, meaning 
that the hailstone will transfer energy to the colder environment 
which then produces a lowering of the temperature. 

During the next cycle, another layer of supercooled water is 
accreted, and the process is repeated, with the radial temperature 
profile resulting from the previous calculations as the initial con- 


dition. 


2.3 Heat Transfers 


When a hailstone grows in a thunderstorm cloud, energy is 
absorbed and released at its surface in several ways. As was discussed 
earlier, latent heat is released by the freezing of the accreted super- 
cooled water droplets; heat is also absorbed in warming the accreted 
deposit, as well as the hailstone itself. Heat is dissipated by 
forced convection and conduction to the ambient air, and latent heat 
is absorbed by evaporation or sublimation of some of the deposit. 

There is also heat generation due to the absorption of the kinetic 
energy of the accreted deposit and due to friction between the moving 


hailstone and the environmental air. Finally, heat can be absorbed 
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or released by radiative transfers. 

In our model the total heat of fusion available is given by 
Cquationm 207) 

ne (OC) (1 - I) 4nR26Ro, (277) 

The energy released by the freezing of the fraction I of 
the deposit has been used up in warming the latter to OCamedhe only 
other heat transfers to be considered in this model are the heat con- 
ducted to the hailstone interior, heat transfers due to conduction 
and forced convection to the environmental air, and the latent heat 
absorbed by evaporation or sublimation of part of the deposit. The 
remaining transfers are assumed to be negligible for reasons which 
we discuss in section 2.4, 

The rates at which the hailstone loses heat because of con- 


duction and forced convection, Q and because of evaporation or 


CGe 
sublimation, Q.., are given byelListe 1963) sandalistectmalwn (1905, 
1967), 
2 8 
Dee = 1.68 ko Re* D (ty - e A (2.8) 
aay 
Qe¢ = aw) 8 Des i, Re* D (ec, - ey) ? C229) 


8 is the roughness factor which is defined as the relationship between 
the heat transfer through a natural (rough) surface and the heat 
transfer through a smooth surface; it is set at 1.5 (List, 1963). Re 


is the Reynolds number. For a sphere falling at the terminal velocity 
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in air it is given by 


Re = a (2810) 


We, and D are respectively the terminal velocity and diameter of the 
sphete monday isathe Kinematic iviscosityook the air.) According: to 


Risteetea ltl 96/)e. (peinikPa)) 
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a 


substituting in (2.8) together with equation (2.4) for The 


Re = 4.19 x 109 RE-5 leauge woo (Que) 


R is the radius/of the sphere (m), i the air temperature in K, and p 
the pressure in kPa. 
k is the thermal conductivity of the air. According to the 


thermodynamic theory for dilute gases (Reif, 1965), 
Keio sae (2.14) 


However, we find that the following formula approximates values of the 
thermal conductivity of the atmosphere, according to the American 
Institute of Physics Handbook (AIPH) 1972, up to the third decimal 


place in the temperature range of interest. 


ROS (eS see ei ee) (2.15) 
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The diffusivity of water vapor in air is calculated according to AIPH 


data also by 


+ 75 
on - 
a ~219 ie, Go), (cm2 5s !) , (2.16) 


although Hall and Pruppacher (1976) suggest 


1.94 ; 
al ee eee | (em? st) , (2.17) 
273 f 


In both equations the pressure is in mb. In SI units equation (2.16) 
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becomes (p in kPa) 


1 
S ay! al ee 
Daz e711 0 Ae D (moes +), (2918) 


C1 9 is the product of several constants and the latent 
b 


heat of sublimation or vaporization. It takes the value 


c = 9.84 x 10% J m2? kPa! and a 8.67 x 10% J m3 kPa 1 


: O Oo ; 
when the air temperature is less than -20 C or at 0C, respectively. 


Between -20°C and Ce varies linearly with temperature from c 
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Upon substitution of all these constants in equations (2.8) 


and (2.9) we obtain 


Ce eatatie a0 (2719) 
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2.4 Other Heat Transfers 


The drag force for a hailstone falling at its terminal 
velocity is equal to its weight, so that the heat generated by 


friction per unit time is given by 
Q Fag OY gs (Qe) 


where m is the hailstone mass and g the gravitational acceleration. 


For a hailstone of density p, = 900 kg m 3 


=) 
2 Sgro (als 
Qep ar 2 9exe 10 eR o (Ww) . (222) 


If the kinetic energy of the accreted deposit was entirely 
transformed into heat, the heat generated during one accretion cycle 


would be 


v2 
K.E. = 4nR2p,6R = (2.23) 


assuming that the cloud droplets have a negligible terminal velocity 


compared to that of the hailstone. Since the energy would be delivered 


during a time t, we can calculate an equivalent rate of release 
2 
Vv 

t | 
2 ae 


= WnR2p, or (2924) 


ke 
In order to estimate the energy absorbed or released because 
of radiative transfers, we investigate the infrared and the visible 


parts separately. For the infrared, assuming that the hailstone and 
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the water vapor in the cloud are perfect blackbodies, the radiative 
flux out of the hailstone would be given by 


Uy oe vue = NP (2.25) 


where o is the Stefan-Boltzmann constant, Ty and Urs the hailstone and 
air temperatures. It has been estimated (Kuo-Nan Liou, 1976) that 
cumulonimbus clouds reflect 80-90 percent of the incident solar radia- 
tion. It seems, therefore, reasonable to assume that 10 percent of 
the solar constant would be an acceptable estimate for the diffuse 
visible light in the heart of a thunderstorm. Secondly, we assume 
that an average hailstone absorbs at most 10 percent of the visible 


light incident on its surface (Hobbs, 1974). Therefore, the net 


radiative flux from the hailstone would be 


ia = heR2(o(Tt - ii 70) %5)ies (2.26) 


where s is the solar constant 1.395 x 1023 Wm 2. 

In order to obtain an idea of the relative magnitude of 
all these heat transfers, calculations were made for a hailstone of 
radius 5.0 x 10 3 m, in the following cloud conditions: Mes = 243K, 
p = 34.5 kPa, Mle = 3.5 x 10 3 kgm 3. The equilibrium temperature 
for these conditions according to the List et al. (1965) model would 


be about 262K. Results are shown in Table 2.1, according to these 


Qep + QQ + We 


- 004 . (2.27) 
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Table 2.1 


Heat transfers for a hailstone 5 x 10 2m in 
radius in the following cloud conditions: 

De 243K pe= Bubs kPa Wier 25a xe Qe kqume. 
Qec is the heat exchange due to conduction and 
convection to the surrounding air; Qpco: heat ex- 
change due to evaporation or sublimation; Qer: 
heat production due to friction; Qa: radiation 
absorbed by the hailstone; Qyp: heat production 
due to the transformation of the deposit kinetic 
energy. 
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Therefore under these conditions and probably under most 
other conditions encountered by a hailstone, it seems to be quite 


reasonable to neglect Ores Q, and Que altogether. 
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CHAPTER 3 
MATHEMATICAL METHODS 


3.1 Two Initial Boundary Value Problems 


The temperature change inside the hailstone is determined 


by the heat diffusion equation 


sr (o,c,T) ~ We(K,¥T) - g = 0 (3.1) 
where p; is the density of the hailstone, C; the specific heat 
capacity, T the absolute temperature, K. the thermal conductivity, 

and g the heat generated per unit volume. According to the assumptions 
made in Chapter 2, PC, 4K; are constant in time and space, g is 

zero (freezing of spongy ice inside the hailstone is not considered 
since the external cloud conditions are kept constant) , and we have 


radial symmetry; the equation therefore reduces to 


oP AL Aen? (332) 
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For further simplification consider the following function 
u = r(T - 273), where u is the product of the radial distance and the 


temperature in degrees Celsius. Upon substitution in equation (3.2) 


we have 


ak (3.4) 


with the boundary condition that u = 0 at r 


ti 
=) 


In practice, the accretion cycle consists of two phases 
which form two distinct initial boundary value problems (IBVP). For 
the first phase, the surface temperature is maintained at 0°C and 
the initial temperature profile f(r) is given. The time duration of 
the phase, te, is not known but has to be determined from the problem. 


Stated formally this IBVP can be written as follows 


du 92u 


=— = | == ; 0 Sir < R, (3.4) 
f ar2 

wi) SS wey ceo), (3.5) 

mip ei Saat) (3.6) 


The solution to this problem is given by Carslaw and 


Jaeger (1959), 
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Since the model considers only conditions where the air 
temperature is below freezing, in order for the hailstone surface 
temperature to be maintained at Onen heat must be generated at the 
surface. We assume that the thickness of the accreted deposit is 
sufficiently small relative to the radius of the sphere for the de- 
posit to be considered coincident with the surface. The source of 
energy is then the latent heat of fusion released on freezing some 
of the liquid water trapped in the deposit. The maximum amount of 
heat available for this is given by equation (2.7). The heat fluxes 


from the hailstone to the air are given by equations (2.8) and (2.9) 


Q(u(R,t)) = Qe + Qe (3.8) 


and the heat flux from the hailstone surface to its interior is given 


by 
, oT du u 
be Die = om pe 
Dyes Mia 5 oye ink ® (Z|, u | : (3.9) 
Therefore, the maximum time duration of the freezing stage 
t pe can be obtained by solving the following equation: 
m 
Uf ( ) 
5 ou UR 
nR*6Ro, (1 - T)L,(0°C) = ha, (ute,e)+ ink R ( 22 a wae) ae 
fe) 
(31.0) 
If we find that 
<2 oT (3.11) 
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where t is the time duration of one full accretion cycle as determined 


by equation (2.1), then 
ti tec (S512) 
If, on the other hand, it happens that 


a) Se (3.13) 


then 


ima, (3.14) 


In this last case the cooling phase does not take place. 


From (3.8) the final temperature profile is 


foe) “ k é 
u(r te) = = E sin(@) exp fsgzee | f Ff (p) sin(= >) dow (3.15) 
n=1 
The time duration of the second phase, the cooling phase is 
given by 


(3.16) 


During the cooling phase, the surface boundary condition is the 
following: the heat transfered from the hailstone to the environment 
must be equal to the heat flowing through its surface. The initial 
temperature profile for this stage is the final temperature profile 


of the freezing stage and is given by equation (3.15). This second 
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IBVP can be written formally as follows: 


oe O<reR (3.4) 
u(O,t) = 0 (3a) 
- nk R? = a = 7 AnK,R nal : HR.) Sige Sh Gh) 
Chas 0) ae a i) ee Le Ce) (3.19) 


Because Qeo depends on the saturation vapor pressure over 
the hailstone surface, which has a non-linear dependence on the tempera- 
ture (exponential in T ! according to the Clausius-Clapeyron equation), 
this second IBVP is not a standard problem with well-known solutions. 
However, Macklin and Payne (1967) did linearize the boundary condition 


(3.18) by setting 


Tiv= ce) : (3370) 


where Te. is the hailstone surface temperature and UL. is the air 
temperature. The proportionality constant he is calculated by 
evaluating Q., and Q.. (they actually used slightly different but 
equivalent expressions) using the deposit equilibrium temperature Ty 


as obtained by the steady-state models and by dividing by the differ- 


ence between this temperature and Te 
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S 3 (352i) 
4k R (T T,) 
With this approximation the solution for the temperature can be 
written as (Carslaw and Jaeger, 1959) 
% ~ka2t R202 + (Rh. - 1)2 
Gegon T+ Fema ie ne ork 
ie O eae s 
n=| Rea 4 Rh, (Rh, iF 
R 
x J (u(o,t¢) - eT.) sina p dp , (3522) 
where the a are Chey nearoo tse of: 
Ra/cot Ra + Rho - 1 = 0. (3523) 


As mentioned earlier in the discussion of the Macklin and 
Payne model in Chapter 1, it is difficult to evaluate the error intro- 
duced by linearizing the boundary condition. Furthermore, it is not 
necessary to look at solutions (3.7) and (3.22) in great detail to 
see that these cannot be of much practical use. Notations like f 
and f are quite elegant, but getting actual numerical evaluations 
may require several approximations or simplifications which would 
have to be repeated for every point in time and space. In this case 
it seems just as well to calculate approximate solutions directly 
from the differential equation (3.4) together with the suitable 


boundary conditions. This is the approach taken in the proposed new 
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model: the partial differential equation (3.4) is approximated by 


means of finite differences. 


3.2 Finite Differences and Grid Mesh 


The simplest finite-difference approximation of the partial 


differential equation is given by the explicit marching form 


Uline iC) surat UCEP EAR, De ce2u (iat) eEruces Ait) 
Se ea eee k te Cae Sky, gama nieaa Es mse2H) 


Thegetoresn | hetheavaluessoh Usatepoints: ter - Ar. Gg toAr are known 
ab time t, tit is possibleito icalculateculr st teat). 

In order to be able to use the finite difference equation 
(3.24), the domain of solution has to be covered by a grid of points 
separated by space intervals At. In the new model the ''size'' of the 
domain of solution is changing constantly from accretion to accretion. 
The are of the spherical hailstone increases with time, and the 
time duration of one cycle t might increase or decrease depending on 
the relative variations of the terminal velocity Vio and the liquid 


water content we (see equation 2.1). Using the conventional notation 


ue = u(mAr,ndt), (325) 


equation (3.24) can be rewritten as 


n n n 
u = Mu oe ue) + (1 - 2h) Lee (3.26) 
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where 


laid a (3.27) 


Fourier component analysis of equation (3.26) yields the 


following stability condition 


(3.28) 


i 
N]|—- 


So once the radial interval Ar has been fixed, there are 
severe limitations on the time interval At. 

In the new model it was decided to use a constant number of 
radial divisions, so that the radial interval Ar would increase from 
accretion to accretion. Preliminary trials were done with the explicit 
scheme (3.26). However, it was found that the method was extremely 
inefficient. Even when the largest Ar possible for acceptable accuracy 
was used, constraint (3.28) was allowing only very small At, result- 
ing in a very large number of time steps for one accretion cycle. 
Furthermore, in the cold parts of the cloud, in the dry growth regime, 


t- is much smaller than i“ with the result that only a few time steps 


i 
were available for the integration of the freezing phase, compared 

with a few thousand for the cooling phase. This was judged unaccept- 
able, since from the heat budget standpoint the freezing part of the 


accretion cycle is just as important as the cooling part, and similar 


accuracy is needed in both. To remedy this problem, it was decided to 


resort to a more efficient implicit method. 
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3.3 The Crank-Nicholson Method 


Another approximation to the partial differential equation 


(3.4) is given by 


Ue eee it) eae (oat) 
At 
ae Pian terete cu (ete te At) erst ae Caan C) (3.29) 


hing 


Here the approximation of the second derivative is expressed at time 
t + At. The Crank-Nicholson method consists in using a linear combina- 


tion of equation (3.24) and (3.29) to obtain the following formula: 


(i= y) uPtl pad = ya ty ut! (1 = y) ult 
= - yAun - (1 - 2ry) u ~ YU : (3.30) 
where y is the weighting factor, and 
Oleeyew lie (323) 


When y = 0 we recover the explicit formulation (3426). 

Upon writing equation (3.30) for all values of m = 0,M, we 
are led to a tridiagonal system of M - | equations in M - | unknowns. 
The reason why there are only M - 1 unknowns at time n + | is because 


ura and ne are prescribed by the boundary conditions. In matrix 


M 


form, the system of equations can be written 
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YAU,» (1 2yA) uy Yuy 7 uy, 
where 
eee (Te -ys (3.33) 
ce Senna) Saber (3.34) 


This tridiagonal system can be solved easily by using the 
Thomas» algorithm. for simplicity let us denote the M >=) 1 "terms; on 


the R.H.S. of system (3.31) by qd: Using the transformation 


( e = z ’ (3.35) 


dé I Ge) (3.36) 
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m i * ’ Wh Pape ci Vy (737) 
oS bie 
m-] 
x an i ad 
Yi EERE m=2,M- 1, (3.38) 
i ele 
m-] 
we obtain, 
n+] us EF 
OM ee] (3.39) 
nh lb? ae sear cant | e , 
u See meg C ae ear mos a2 ee (3.40) 


This scheme is stable for all A (Forsythe and Wasow, 1960). 


3.4 An Iterative Boundary Condition 


During the cooling stage of the accretion cycle, the amount 
of heat transferred to the atmosphere must be equal to the amount of 
heat flowing through the hailstone surface. This is the external 
boundary condition of the second IBVP, 


du Meu Rec) . 
: bak r (22 Rae ) Dep tls © (3.18) 


Writing this in finite difference form, 
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n+] _ n+] a 


ey uy = Ar/R)eeeRATO (3.42) 


From equation (2.9), Q., is directly proportional to the saturation 
Vapor pressure over the hailstone surface, which depends on the surface 
temperature in a highly non-linear fashion. Equation (3.42) cannot be 
solved analytically in terms of ve to permit the use of the Thomas 
algorithm. Instead, an iterative method has to be used. Furthermore, 
it was found that a second order approximation of the surface derivative 
a Rp gave better results; nerene this is because the finite differ- 


ence approximation of the heat equation itself is of second order in r. 


Hence, we write: 


n+] n+] n+] 
3u See ae a 
oe 2 = OED SSN ee RON? ae O(Ar2) . (3.43) 
h 2Ar 


Upon substitution of this second order formulation in equation (3741) 


we get 
n+] n+] n+] n+] 
3u sll O Vs pe edb l ae ». 
Lee eile er ayy? (3.44) 
2A R 


Using the Thomas algorithm to get Uno and Usps substituting in 


(3.44) and rearranging, we get 
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A = YAU, _» (1 2 hy) Uy] - YAUy : (3.46) 
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Since Qo” is actually a function of u ,» eguation (3.45) contains only 


M 


n+] 
one unknown u 


ft and it can be solved iteratively. In the model this 


is done by the method of false position (Sokoinikoff and Redheffer, 


1966). Denoting the L.H.S. of (3.45) by F(x), the root of 


EWE SW) (3.47) 


n+] 
M 


as" , 
would be the proper value for Us lb Since the upper limit for u 


is 0 (corresponding to the maximum temperature allowed in the model, 


0°C), it can easily be shown that the following recurrence formula 


x ,F (0) 


Xx 
n+] SCOre FOr (3.48) 


n 


will converge to the root of (3.47). As a first guess for x, the 
value corresponding to a surface temperature equal to the air tempera~ 


ture is used since this is the lowest possible value of the surface 


temperature. 
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3.5 Stability with Respect to Heat Transfer 


The duration of the freezing phase is determined by the 
amount of latent heat available from the freezing of the deposit. 
The time necessary for complete freezing is found by solving equation 


(S210) 


mf 


t 
WmR?6Ro (1 - T)L-(0°C) = f ay(u(R,t) + bane ( Se, wet) ) ae 
Oo 


Since we are using finite-difference methods, the R.H.S. of 


equation (3.10) becomes a summation: 


f 
Me 


(Note that here N is the number that will make (3.49) equal to the 


i a) es 


ul" bs) u” fe 
Q, (ur) + 4nk.R ea - = At . (3.49) 
ty i R 


n 


L.H.S. of (3.10).) Unfortunately expression (3.49) proved quite diffi- 
cult to handle as it was found to be very sensitive to the magnitude 
of Ar and At. This is most likely due to the behaviour of the sur- 
face temperature gradient at the beginning of the phase. Initially, 
the surface temperature is raised suddenly at Once this leads to an 
extremely large surface temperature gradient which decreases very 
rapidly afterwards. If an average of the initial and final tempera- 
tures for one time-step is used in an effort to get a better approx- 
imation to the average gradient, the estimated freezing time (NAt) 
changes unpredictably when Ar and At are made smaller independently. 


This is basically due to the fact that the gradient itself decreases 


- ha 
a ve Map lytek Hi ery ea ae By al ll 


Rs an! oe see? otal, 
Mi Me Hel GC Ay) 8 pea ‘Ves ig? oeeaiee a 


3 « 
7 . 
~ 7 oy + = ‘ 
, i 2 a" ‘ “ 
: i j 
4 as Me’ es a . @,] : y De ine «6 1) 
oo 4 ae oe 7 ye Biss 
Toe @ 
~~ S. z= 
: y Wee 4 ‘ - 
5 a a Fe 
° $F My r q - 
L on s 4% 
# 
é 
. r 7 e bs n ) - - 
iq ' 6 .@ 
14 ; , : : 
* é ws - ead 
_ a & “ 
— on * cia pi . - yt! & ~ - = eal al 
ms ¥ - i : ¢ , ae 
—. a my ay @ os ‘i if 
f ¢ Wy ‘ 
Hk) @ Th eit erie mr 1608"; ‘ plalis epg ay 4g af 
a 2 2 1 “ ae eee 
lb ained Ai 1) a | 
to + - 
" Si AR SiS) eprige’ as /65 hc ds eps Thee el @ Pi 
tall 2 & a i S ‘a ba ia s 4 
‘ 
avi Peter’ yapoan ba | wu feat af? et d yy iets 
~ ie ~ > i PS Sane = 


hd) Ghia Ah: app Pee tei .« ATE “ 
es ee a = ~ mw Ee & fino Pn = 


view. teiel..c Gpiee: peBENY A POG OEE NY Gy, 8) 90: 
ee de ee “APs AY. x hes 
ioe ring} ia pe? i a’ ee ai ud Sols te 
wit ras e ie PES : : bade i ae oe 


rey oe ar ba Het ee 


43 


exponentially with time (it can be easily seen that differentiation 
of (3.7) with respect to r does not affect the time dependence); 
therefore a simple average of the initial and final gradients, if 
those were exactly known, would always give an overestimate of the 
average gradient. On the other hand, in this case, the first order 
finite difference approximation always underestimates the surface 
gradient: so depending on the relative magnitude of Ar and At the 
internal heat fluxes would be underestimated or overestimated. A 
similar situation occurs when only the initial temperature for one 
time step is used. If the final temperature alone is used, the 
internal heat fluxes are always underestimated; but in order to attain 
a sufficient accuracy, Ar and At have to be so small that the ad- 
vantages of using the Crank~Nicholson method are lost. In the search 
for a better way to estimate the internal heat fluxes, it was noticed 
that the temperatures themselves did not seem to be significantly 
affected by changes in At and Ar. It was therefore decided to 
calculate how much the total heat content of the hailstone changed 
from one time-step to another. That is, at each time step, the 

total heat content is calculated and compared to the initial heat 
content to see how much heat has been conducted inside. This tech- 
nique works quite nicely and although it requires more computation 

at each time iteration, this is more than compensated for by the 
fewer iterations needed with the Crank-Nicholson scheme. The in- 
creased accuracy of this method might also be explained by the fact 
that the total heat content is essentially calculated by adding the 
temperature at all the grid points, and its relative error remains 


of the same magnitude as that of the temperatures, whereas the 
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relative error of the surface gradient, which is essentially a differ- 
ence of the temperatures can become intolerably high. (Note that 

even if the gradients might be very large due to very small Ar, the 
difference between two adjoining temperatures can be quite small.) 

In the solution of the second IBVP of the second phase, when 
the change of the total heat content of the hailstone was compared to 
the heat that should have been transferred to the environment accord- 
ing to the surface temperature calculated from the boundary condition 
(3.43), a similar discrepancy occurred. Independent changes of Ar 
and At resulted in surface heat transfers, as calculated from (2.8) 
and (2.9), that were sometimes smaller, and sometimes larger than the 
internal heat changes, again in a rather unpredictable manner. JHow- 
ever this time there seemed to be no possibility of avoiding the use 
of the surface temperature gradient. Although this problem might look 
similar to the one encountered in the first IBVP, it is different in 
the sense that we are not trying to estimate heat transfers from sur- 
face gradients during the span of one time-step. The surface tempera- 
ture is found by equating the heat flux due to the gradient of 
temperature at the surface, to the fluxes to the atmosphere due to 
the surface temperature, at one particular instant. Although there 
may be some error in the approximation of the surface temperature 
gradient, the main cause of this problem is likely due to the method 
used to integrate the heat diffusion equation itself. 

Indeed, the basis of the Crank-Nicholson scheme is that the 
rate of change of the temperature at one point during the time interval 


At is set equal to the average of the finite difference approximation 
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of the laplacian of the temperatures expressed at time t and at time 
t + At; this should give a better representation of the average rate 
of change of the temperature during the time interval At, than the 
simple explicit formulation (3.23). However, from the theoretical 
solutions (3.7) and (3.21), we see that the temperature has an 
exponential ''decrease'' with time. Therefore as the length of a time 
step increases, a correspondingly greater weight should be given to 
the laplacian expressed at time t + At, in order to have a more 
accurate representation of the average rate of change of the tempera- 
ture during the interval At. The original Crank-Nicholson method 
gives equal weight to the two laplacians. In order to remedy this 


problem the weighting factor y was set equal to 


(es seat ee (3.50) 


where A is given by equation (3.26). Therefore as At increases, y 
decreases and the weight is shifted towards the laplacian at time 

t + At. When At decreases, A becomes smaller, and y tends towards 
1/2, that is, equal weight is given to both laplacians. This method 
totally eradicated the wild fluctuations between the surface heat 
transfers and the internal heat changes, altnough some consistent 
error persisted. This error could be reduced at will by reducing 

Ar and/or At. Since no improvement was obtained by expressing the 
surface gradient by an approximation of higher order than second 
order, it was decided that the remaining errors were simply truncation 


errors of the finite-difference approximation. Maybe an improvement 
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-could be attained by devising a more clever y, or by simply going to 

a more accurate scheme for the integration of the heat diffusion 
equation. However the present method is judged accurate enough for 
our purposes. Furthermore, according to Forsythe and Wasow (1970, 

p. 124) this method is stable and convergent, no matter how At and 

Ar tend to zero, when the weight given to the laplacian at time t + At 
is equal to or greater than 1/2. This was amply verified in the many 
trials done to find a suitable combination of the number of time and 


Space divisions. 


Ne ee ae 
- soe 
a ise ek 
SORT! <ehe U: Rae a8 “ETE, = 


ar od repeal Ji mere ing ea 
ne ee oe! '@i Farts ng alia 


Vode 004 ed ‘ ry wae ” vere ie woe 


Se). 60s >7 fii tse ete ons = ae 


4 - 
= . 
#i 
- “ 
ry a t 
a 
= S - + es ” 
- ~ i Gt . ‘ 
. - % 
; : =f, 7 a ie ie 
: g 
a 

x = . “i 
H me 
* 

> ee . + . y 
‘ 
- a - i“ 
¢ 
— = = 7. Gets 2 se 
a 
- * is = “<1 7 - Lr a 
Pe AT ee We - 
rl : 
~ P - “ ~ 
~ € * . ax Fe ou A = 
- aed + 
e ¢ . oer 


CHAPTER 4 


A FEW PROGRAMMING CONSIDERATIONS 


4.1 Brief Description of Program 


As previously mentioned, the model investigates the growth 
of a spherical hailstone under constant cloud conditions. The program 
is divided into three parts. The main program calculates the terminal 
velocity and the time duration of one cycle according to formulas 
(2.1) and (2.3). It also performs various other tasks like the 
calculation of the resulting ice fraction, the re-initialization of 
the grid after one cycle and the output of results. The solution to 
the heat diffusion equation is handled by two subroutines, FREEZ and 
COOL, depending on whether the hailstone is in the freezing stage or 
the cooling stage of the cycle. The mathematical methods used in 
these routines have been described in the previous chapter. A copy 
of the program is provided in Appendix 2. Calculations for LO 
accretion cycles required about ten seconds of computing time on the 


Amdahl 470 at the University of Alberta. 
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4.2 Re-initialization After One Cycle 


In this model, the number of grid points along a radius 
remains constant from cycle to cycle, and since the radius is becoming 
larger, the distance between two grid points is also increasing. 
Therefore from one cycle to another, the grid points do not correspond 
exactly. At the beginning of a cycle, immediately after the accretion 
of a new layer, the temperature is not known on the expanded grid, and 
it must be interpolated from the final temperature profile of the 
previous accretion. Let e be the displacement of a new grid point m, 
from its position before the accretion; then the temperature at this 


point Bue can be approximated by: 


Uo = Ut |, (4.1) 


where u is the temperature at point M in the previous grid, and Ar 
is the previous grid length. The displacement ¢« is given by: 
6R 
€ = ore 6 (4.2) 
where M is the total number of divisions, m the grid point number, and 
SR the thickness of the accreted layer. At the surface point the 
temperature is set to Once according to the model assumptions. 
Some difficulties may arise however when the deposit thick- 
ness 6R is of the same order of magnitude as one space intenvalenreaenld= 


ure 4.1 illustrates schematically what can happen when 6R is slightly 


greater than Ar. In this figure we see that the value at grid point 
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Figure 4.1 Problem of grid re-initialization when §5R the deposit 
thickness is larger than Ar the grid point spacing. 
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7 (or 8) in the new grid would be better approximated by using esti- 
mates of the slope at grid points 8 (or 9) in the old grid, instead 
of at points 7 (or 8) as per equation (4.1). Nevertheless when the 
situation arose, total heat contents before and after the re- 
initialization were found to be within less than .05 percent of each 
other, indicating that procedure (4.1) was accurate enough for our 
purposes, 

On the other hand, grid point 9 is now inside the accreted 
layer. In the model, if this situation were to occur, this grid 
point would also be initialized to 0°C (indeed it was decided that 
whenever a grid point would be within .2 Ar from the edge of the 
accreted layer, in the new grid, this grid point would also be 
initialized to 0°c); however during the integration of the heat 
diffusion equation for the freezing stage only the end point would 


be kept at 0°C as a boundary condition. 


4.3 Ice Fraction 


When a hailstone enters the wet growth regime, the surface 
deposit freezes incompletely, because the surface heat fluxes are 
insuftticitent to get rid of all of the latent heat. Siiseresults in 


a new ice fraction L' which is calculated in the following manner: 


Tien ee ae) eee / Lmao | (4.3) 


where I is the initial ice fraction of the deposit as given by ORL 


Hoot is the total heat that has been conducted inside the hailstone 
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and transferred to the environment, and eee is the maximum heat 
which could be available if the deposit were to freeze entirely. 

In this regime, the surface temperature remains at 0-ce 
and the hailstone interior gradually warms up towards the melting 
point. After a certain time the internal temperature becomes uniform 
enough that it is meaningless to solve the heat diffusion equation. 
At this point the internal temperature gradients are so weak that 
the amount of heat conducted inside the hailstone becomes negligibly 
small compared to the other heat fluxes. The ice fraction is then 
basically determined by the external heat transfers alone and the 
model becomes essentially similar to the continuous models discussed 
in Chapter 1. It was decided (more or less arbitrarily) that this 
point is reached when the hailstone average temperature is equa! to 
or greater than -.001°C. Since the hailstone cannot become warmer 


nN 


than 0°C this would correspond to a very weak gradient indeed. 


h.4 Heat Content 


As was discussed in Chapter 3, the easiest way to calculate 
how much heat has been conducted inside the hailstone is to evaluate 
the change in its total heat content. In our case the total quantity 


of heat H. contained in the stone is given by 


R 
Hews f c. T(r) harp dr (4.4) 
Oo 


if we express T(r) in degrees Celsius and make use of the function 


u = rT, we can define 
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R 
H. = AT .C. j ru(r)dr . (4.5) 
fo) 
Therefore the heat content of a hailstone at 0°C is ZEhO wean Celts 
negative for lower temperatures. In view of the initial assumptions, 
p; and c, are constant throughout the hailstone. The temperature is 
known at every grid point and integral (4.5) can be accurately approx- 


imated by the trapezoidal rule 


R = 
f u(r)rdr = (Ar)? © mu + vu (4.6) 
oO = 


An interesting paradox arises in the calculation of the 
initial heat content, at the beginning of the freezing phase. Since 
in fact, the model accretes a deposit which is already at 0°¢ (the 
heat of fusion released by the initial freezing having been used to 
warm the deposit), the total heat content of the hailstone is not 
changed, according to (4.5). However since the total mass of the 
hailstone has increased, the mean temperature has also increased and 
the hailstone is in fact warmer by the amount of heat that has been 
released during the initial freezing. Nevertheless this effect will 


have to be taken into consideration in the calculation of heat budgets. 


4.5 Saturation Vapor Pressure 


tn order to calculate the heat transfers due to sublimation 
or evaporation of some of the deposit, it is necessary to calculate 


the vapor pressure at the hailstone surface and in the environment. 
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The model assumes that in the cloud the air is saturated with respect 
to a plane water surface at the environment temperature. This vapor 
pressure is calculated by Richard's method as formulated by Wigley 

(1974). The vapor pressure at the hailstone surface is assumed to be 
the saturation vapor pressure over a plane ice surface at the surface 


temperature, which is calculated by Lowe's polynomial approximation 


Go7n) 


4.6 Thermodynamic Constants 


The value of the thermal diffusivity Kk. used in the inte- 


gration of the heat diffusion equation is 1.09 x 10 2 cm TCalaset 


(since the program uses c.g.s. units) as proposed by Picca (1964); 


Cum 


in S!_ units this is 1.09 x 10 +m <C ae 


Ss The value of the thermal 


conductivity K, selected for the hailstone is 2.09 x 10 * Wm } tee 
(.00500 c.g.s.) also as given by Picca. Measurements (Hobbs, 1974) 


indicate that it could vary between 2.5 x 10 * wm! %?! and 


2.0x 10 *Wm! °c ! over the temperature range -30°C TOMO Chany 
value of 1.46 x 10% Wm! °c! (.0035 c.g.s.) was used by Goyer et 
al. (1969); however they claim that variations "over the physically 
reasonable range showed that the results were quite insensitive to 
this parameter". 

Even though the heat capacity of ice can vary as much as 
17 percent over the range of 0 - CORON we are constrained by the 
assumptions of the model to keep it as a constant. In order to.be 
consistent with equation (3.3) and the values selected above, the 


heat capacity of our theoretical hailstone was set at 2.09 J kg ! 
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(7500nca leon) 5 503) call qu seis ethe specific heat capacity for ice 


at 0°C given by the Smithsonian Meteorological Tables (List, 1958). 
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CHAPTER 5 


RESULTS AND DISCUSSION 


Dplmecenera Desc ption 


In this chapter, the model is used to investigate the growth 
of a hailstone in constant cloud conditions. Forty layers of 25 um 
each are added to a spherical hailstone of radius 4.5 mm, resulting in 
a sphere of radius 5.5 mm. The details of these growth segments are 
calculated for six different sets of cloud conditions: for the first 
five sets, the stone is in the dry growth regime, and for the last 
one it enters the wet growth regime. These conditions are summarized 
inn ab les5. |. 

T., the initial constant temperature of the hailstone, was 
selected to correspond approximately to the equilibrium temperature 
Us as obtained from the continuous models for a hailstone of 5.0 mm 
radius. It can be seen from Figure 5.1] that the hailstone exhibits 
a fairly similar behaviour in all of the first five sets of cloud 
conditions. The final surface temperature for each cycle drops 
sharply after the first few accretions, then rises gently afterwards. 


Around accretion number 15, the hailstone seems to have reached some 
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Table 5.1 


Environmental conditions for the investigated growth 
segments. P is the air pressure, T. the air tempera- 
ture and we the liquid water content. These conditions 
are consistent with the cloud model used by List et al. 
(1965). Te wise the ini tialeconstang temperatures oietne 
hailstone and T, is the equilibrium temperature as 
obtained from List et al. (1965, 1967) model for a 
hailstone of 5.0 mm radius. The number in parenthesis 
ise the ice fractionr 


P(kPa) = T(°C) sw, (1073 kg m°8) mee) Cp) 
34.5 30) 3.35 171-0 =11.0 
38.0 25 Shee) -8.0 -8.0 
WO V0) 3e5 -6.0 -h.0 
16.0 at) Bau 25 =250 
51.3 -9 2.0 “ioe “1.5 
55.5 -5 2.0 -1.0 0.0.75) 
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Figure 5.1 
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Number of Accretions 


Final surface temperature of the hailstone after N 
accretions for different cloud conditions: an by 

c, d, e, (see Table Couly) m eNote that continuous lines 
were drawn for ease of representation, each curve 

is obtained by joining a set of discrete.points. 


Bu 


. 


2 = eo 
Oo ee = 

ie: ‘Mere s 

e 
-, wl : 4 > a Se i 

. 

. t. er 

a 

7 

1 
Ld ou bd y= 


= e = 7 


: F caine ere 


~~ . fs * , 


+ 


58 


kind of equilibrium with the environment, and the slow rise of the 
final surface temperature can be attributed to the slowly increasing 
radius. 

The slightly different aspect of curve c can be explained 
by the fact that the initial temperature T is 2°C cooler than the 
equilibrium temperature. For this reason the temperature starts to 
rise ''sooner'' than in the other cases, and ends up higher than the 
initial value after the full 40 accretions, also in contrast to the 
other cases. 

For condition f, the hailstone enters the wet growth regime 
after the first accretion (Figure 5.2). The ice fraction drops sharply 
for a few accretions, again some sort of equilibrium after accretion 
number 10. Here also we did not select the equilibrium temperature 
as the initial condition in order to see how the equilibrium would 
be achieved, if this were the case. 

At this point it should be made clear that the time scale 
foretheurealizatlon- or al bithe: cunves inn! Guree>. mandan| gure. ons 
is different for each case. Because the hailstone is in different 
environmental conditions, its fall speed and its accretion rate are 
different. Table 5.2 gives a comparison of the elapsed time after 
the accretion of an equal number of layers for each case investigated. 

Note that the time interval for the accretion of one layer 
decreases as the number increases; this is because a large hailstone 
has a faster rate of accretion than a smaller one, under the same 
environmental conditions. In Table 5.3, the final surface temperature 


for the twentieth accretion cycle is compared to the equilibrium 
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of 25 um layers, for cloud conditions f. 
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Table 5.3 


Comparison of the final surface temperature iiee after 
the 20th accretion cycle with the equilibrium tempera- 
ture T, for a 5.0 mm hailstone under different cloud 
conditions. The number in parenthesis is the ice 


fraction. 
Cloud 
Pena cione 2 : C q e f 
TS ¢ -12.81 -9.62 -5.85 -3.97 -2.54 O70) 
4 -11.0 -8.0 -4.0 -3.0 -2.0 0(.75) 
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temperature as obtained from the Listogram (Appendix 1) for each set 
of environmental conditions. 

It is easily seen that for most cases Toe is significantly 
lower than the equilibrium temperature. The difference however de- 
creases at higher temperatures. For conditions f both models result 
in the wet growth regime, and the ice fractions obtained are in 
reasonable agreement, considering the accuracy with which it can be 


read from the Listogram (Appendix 1). 


5.2 A More Representative Temperature 


The behaviour of the final surface temperature for each cycle 
does not permit one to draw any conclusions about the processes simu- 
lated by the model since this temperature is not representative of 
the hailstone during one accretion cycle. According tolFigure: 5.3 
and Figure 5.4, there is quite a bit of variation near the surface, 
especially for coid conditions. 

Figure 5.3 indicates that after the freezing phase, heat 
has penetrated further towards the center of the hailstone in warm 
conditions than in cold conditions. This is because it takes longer 
for the deposit to freeze in warm conditions. Since the surface is 
maintained at 0°C for a longer period, there is more time for heat 
to diffuse inside. After the cooling phase, at the end of the 
accretion cycle, the temperature near the hailstone center is slightly 
warmer than at the beginning by typically a few tenths of a degree. 


However the surface temperature is significantly lower, as was 
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Figure 5.3 Internal temperature profiles for the first accretion 
cycle, after the freezing phase (those converging to 
90°C) and after the cooling phase for cloud conditions 
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observed in the previous section. It is interesting to note how 
similar in shape are each of the final temperature curves. 

Figure 5.4 shows the same temperature profiles for the 
twentieth cycle. Here the temperature near the center remains vir- 
tually unchanged. Furthermore the depth of penetration (that is, the 
distance from the surface after which there is no significant tempera- 
ture change) decreases for colder conditions. However it is still 
large enough to necessitate the consideration of spherical effects 
in the integration of the heat diffusion equation, which effects were 
neglected by Macklin and Payne (1967). Again the final temperature 
profiles are very similar in shape. 

Figure 5.5 confirms our previous findings about large 
temperature fluctuations near the hailstone surface. We see that 
for conditions a, the hailstone spends about 65 percent of the time 
at a temperature lower than the equilibrium temperature is obtained 
by the continuous models. As the cloud conditions become warmer, the 
hailstone spends a larger proportion of its time at a temperature 
above the corresponding lg. 

At this stage the need for a temperature that is more 
representative of what happens during one cycle is quite evident. 

For this purpose it was decided to calculate a time average of the 
temperature at each grid point for the duration of a cycle. Figure 
5.6 gives the profiles of such average temperatures for the first 
accretion. It is quite noticeable that all the curves are fairly 
uniform, although exhibiting generally a slight fall towards the 


center of the hailstone. Aiso all curves with the exception of 
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Figure 5.6 Time averaged temperature profiles during the first 
accretion, for ali of the dry growth conditions. 
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Curve c (because of the lower ein conditions) have a slight max- 
imum which is close to the surface, but away from it, indicating that 
the surface heat flux was on the average out of the hailstone. 

By the time the hailstone reaches the twentieth accretion 
cycle, its average internal temperature is very uniform. From Table 
5.4 we see that in all cases there is a slight temperature drop from 
the surface to the center, about one-tenth of a degree, indicating 
that the hailstone is warming slowly on the average because of its 
increasing radius. However, essentially, the hailstone has a constant 
average temperature which is indistinguishable from the equilibrium 
temperature obtained by the continuous model. Therefore it seems as 
if the hailstone adjusts itself to a new thermodynamic equilibrium 
in such a manner that its average temperature is the same as the 


temperature predicted by the continuous models. 


5.3 The Mechanism of Equilibrium 


lt was mentioned in section 5.1 that after the first 
accretion, the final surface temperature is significantty lower than 
the initial temperature; also this effect is more pronounced for 
colder conditions (see Figure 5.1). Similar results were obtained by 
Macklin and Payne (1967), although they seemed somewhat surprised 
by these findings. They thought that the cooling phase would terminate 
when the surface temperature returned to the initial temperature, and 
that this would happen when all the heat conducted jaternally during 
the freezing phase had been transferred back to the environment. 


But this is impossible and the reason is very simple indeed. Since 
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Table 5.4 


Average temperature (20) at different points along the 
hailstone radius, for the twentieth accretion cycle, 
under different cloud conditions. Peis thes coprespond= 
ing equilibrium temperature from List et al. (1965) 
model. 
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by hypothesis, the surface temperature is raised to 0°C during the 
freezing phase, this means that part of the hailstone interior will 
warm up and the internal temperature will rise (see Figure 5.3). When 
the hailstone enters the cooling phase, its surface temperature falls 
faster and becomes lower than the temperature at some point inside 
(otherwise there would be no heat flow out of the hailstone).. Asa 
result, when the hailstone surface returns to the initial temperature, 
at some point inside, the temperature has to be warmer than the 
initial value. It is therefore physically impossible for all the heat 
that has been injected into the hailstone during the freezing phase 

to have been transferred back to the environment at that point. 

We find that the surface returns to the initial temperature 
fairly early during the cycle (Figure 5.5 is a good indicator al- 
though not exactly representative of the first cycle) and that the 
hailstone ends up colder than initiaily after the cycle is fully com- 
pleted. This is because the amount of heat conducted into the hail- 
stone during the freezing phase is significantly lower than the amount 
of heat transferred back to the environment during the cooling phase. 
Intuitively, however, we know that if the external conditions are 
kept constant, the hailstone cannot keep cooling forever, accretion 
after accretion. There must come a stage where the heat transfers 
will adjust themselves to try to balance on the whole. 

In Table 5.5 we see how this process takes place. During 
the first few accretions there is a net cooling of the hailstone which 
causes a lower final surface temperature; this in turn causes sharper 


internal gradients at the beginning of the next accretion resulting in 
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Table 5.5 


Energy budgets of each accretion cycle for cloud condition 
a. Hj is the heat conducted to the hailstone during the 
Freezing stage; Hy is the heat transferred to the environ- 
ment during the cooling stage; Hq is heat released by 

the cooling of the deposit (see text); W is the net warm- 
ing of the hailstone. 
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a stronger internal heat flow. Therefore it becomes apparent that 
the hailstone attempts to adjust its temperature cycle in order to 
reach a final temperature profile such that the heat which flows in- 
side during the next freezing phase will be balanced by the heat that 
flows out during the subsequent cooling phase. Unfortunately the pro- 
cess of achieving equilibrium is further complicated by the fact that 
the hailstone is growing and that the accretion rate (and therefore the 
rate of heat release), increases with the hailstone radius faster than 
the cooling rate (see Chapter 2). This is why there is a point of 
turnaround near accretion number 20 and why the hailstone warms up 
slowly thereafter. 

We will now explain in greater detail how the budgets of 
Table 5.5 were calculated. H. is the amount of heat that was con- 
ducted into the hailstone during the freezing phase; it is equal to 
the latent heat released by the freezing of the remaining water in the 
deposit after the initial freezing, and is given by equation (297)R 
H is the amount of heat that has been transferred from the hailstone 
to the environment. This can be calculated by using equations (2.19) 
and (2.20) in conjunction with the surface temperature at every time- 
step, or by evaluating the difference in the hailstone heat content 
at the beginning and at the end of the cooling phase. Both methods 
were used and it was found that values obtained with the first method 
were consistently higher than the values obtained by the second method 
by a little less than 2 percent. However, as mentioned in the dis- 
cussion of section 3.5, this was judged to be an acceptable error. 


A slight difficulty arises because of the procedure used to 
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calculate the hailstone heat content. The total heat content of the 
hailstone remains unchanged when the deposit is instantaneously 
accreted on it, although the hailstone is in fact warmer (because of 
a higher average temperature as discussed in section 4.4). Since the 
deposit is maintained at Ome throughout the freezing phase, none of 
its heat could diffuse into the stone interior, and only the latent 
heat released has affected the total heat content. However, during 
the cooling phase heat is removed from the deposit also, and this is 
debited against the total heat content of the hailstone. Therefore 
in order to obtain the net warming, the heat released by the cooling 
of the deposit, Hy in Table 5.5, must either be added to the heat 
conducted to the hailstone (H.), or subtracted from Ho 

Table 5.6 gives the net hailstone warming for one cycle at 
a selected number of accretions for all other dry growth conditions. 
Ciearly the hailstone exhibits a similar behaviour in each case (with 
the exception of case c) and by the twentieth accretion an equilibrium 


fsvattained. 


5.4 An Hypothesis 


Let us return to Figure 5.5 and examine one curve in particu~ 
iar, say curve b. This curve registers the surface temperature of the 
hailstone for the twentieth accretion cycle. However, instead of 
assuming that the deposit has been accreted all at once, let us con- 
sider that it is being formed of individual cloud droplets arriving on 
the surface at short regular time intervals. A droplet that had just 


arrived on the surface would have just entered the freezing phase and 
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Table 5.6 


Net hailstone warming in Joules, for one accretion 
after a selected number of accretions. These are for 
cloud conditions b, c, d, and e (see Table 5.1). 
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would have a temperature of o°C; but a droplet that had arrived one 
full period t earlier would presumably be at the end of the cooling 
phase and would have reached the final temperature. Similarly a drop- 
Tet accreted half a period ago would have cooled to a temperature of 
about -8.5°C. Since on the average, droplets should cover the entire 
hailstone surface within one period t, one could make the hypothesis 
that curve b represents the areal distribution of temperature for a 
hailstone of 5 mm radius, at equilibrium, growing in environmental 
conditions b (see Table 5.1). Of course such reasoning could also be 
applied to curves a, c, d, and e. 

Unfortunately, there are several objections against this 
hypothesis. The first one, and undoubtedly the most important, is 
that our calculations assume radial heat transfers. If droplets are 
to be considered individually, lateral heat transfers would have to 
be taken into account: so one cannot just simply assume that individual 
droplets would have a temperature cycle identical to that of a uniform 
deposit. A second objection lies with the statistical nature of the 
accretion process. It is very unlikely that the hailstone will be 
covered uniformly by the droplets; more likely droplets will overlap 
in some portions of the surface while they will be late in covering 
some other regions. (And we won't even think about possible splash- 
ing or bouncing of the droplets!) 

it would be very difficult to assess quantitatively the 
total effect of these processes. Considering the freezing phase, for 
example, it would appear that the lateral heat transfers might shorten 


the freezing time. On the other hand, the presence of several droplets 
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at different temperatures over the hailstone should have a tendency 

to smooth out the internal temperature gradients; this could have the 
Opposite effect of reducing the freezing rate. Alternatively, during 
the cooling phase, a droplet could be affected differently depending 
on whether it happens to be in an area where the majority of the drop- 
lets are at a colder or warmer stage of their cycle. 

Keeping those objections in mind, let us nevertheless go 
ahead with the hypothesis that Figure 5.5 gives the areal temperature 
distribution of a growing hailstone in equilibrium at different cloud 
conditions. Several interesting characteristics of the hailstone 
thermodynamic state = be estimated from this hypothetical tempera- 
ture distribution. The most obvious is of course the mean surface 
temperature; this in fact corresponds to the time average of the 
surface temperature already discussed in section 5.2 and found to be 
essentially equal (within a few tenths of a degree) to the equilibrium 
temperature obtained by the continuous models. 

A second important parameter of the temperature distribution 
is the standard deviation. It gives a better idea of the spread of 
the surface temperatures about the mean, which is not quite obvious 
just by looking at Figure 5.5. Such standard deviations were cal- 
culated for the five cloud conditions leading to dry growth for a 
S mm hailstone and are plotted in Figure 5.7. When the hailstone 
reaches the wet growth regime, the model yields a uniform temperature 
of Ogce with a zero standard deviation. This is why the curve goes 
through the origin. An interesting feature of this graph is the 


weak dependence of the standard deviation on the mean surface tempera- 
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Standard Deviation (°C) 
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Mean Surface Temperature (°C) 


Standard deviation of the surface temperature distribution 
as a function of the mean temperature of the hailstone 
surface for a 5 mm radius hailstone. 
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ture. For temperatures colder than -2°C, the standard deviation is 


By calculating the ratio of the freezing time to the total 
duration of one cycle, an estimate of the proportion of the surface 
which is in the wet growth regime can be obtained. This has been 
plotted in Figure 5.8. It can be easily seen that the ''fractional 
wet growth area'' has a very strong dependence on the surface tempera- 
ture. These results would also support the suggestion that the 
transition from the dry to wet growth regimes is not discontinuous 
as implied by the continuous models, but occurs relatively smoothly 


as the mean surface temperature rises towards the melting point. 


5.5 Effect of Deposit Thickness 


ln order to investigate possible effects of the thickness 
of the accreted deposit, the calculations were done again using deposits 
of thickness 5, 10 and 50 um for cloud conditions a and d. !t was 
found that in each case the hailstone undergoes a temperature cycle 
fairly similar to the one described in the previous sections. Further- 
more, according to Table 5.7, the hailstone soon reaches a state of 
equilibrium in which the average temperature rises gradually due to 
its slowly increasing radius. Also, at equilibrium, the average 
temperature is essentially equal for all layers, the maximum differ- 
ence being less than 1°c from the largest to the smallest layer 
under condition a. However, according to Table 5.8 there is a sig- 
nificant dependence of the final surface temperature after one 


accretion cycle upon the deposit thickness. 
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Fraction of Surface 


Figure 5.8 
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Mean Surface Temperature (°C) 


Fraction of the surface in the wet growth regime as a 
function of the mean surface temperature. This is 
obtained by calculating the ratio of the freezing time 
to the total time of one accretion cycle for a hailstone 


of 5 mm radius. 
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Table 5.7 


Time averaged surface temperature over one accretion 
cycle, at selected hailstone sizes, for different deposit 
thicknesses, under cloud conditions a and b. The start- 
ing radius is 4.5 mm for the 25 and 50 um layers and 

4.8 mm for the 5 and 10 um layers. 


OR 5 um 10 um 25 um 50 um 
R (mm) (a) 
4.55 = = -10.93 -9.87 
4.80 -10.63 -10.40 -11.01 -10.96 
4.90 -10.34 -10.61 -10.94 -10.96 
5.00 -10.08 -10.47 -10.82 -10.90 
5.10 -9.89 -10.32 -10.70 -10.80 
5720 AS PS} -10.19 = 10m op -10.68 
(b) 
4.55 . - -2.07 -1.64 
4.80 -2.2] -1.94 -2.28 =i deh | 
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Table 5.8 


Final surface temperatures after one accretion cycle 
at selected hailstone sizes, for different deposit 
thicknesses OR, under cloud conditions a and b. The 
starting radius is 4.5 mm for the 25 and 50 um layers, 
and 4.80 mm for the 5 and 10 um layers. 
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Table 5.9 compares the accretion of deposits of different 
thicknesses for an equivalent radial increment, starting with the 
same initial conditions. When the hailstone reaches 4.60 mm, just 
about the same amount of heat has been conducted to the interior, re- 
gardless of whether the accretion consists of two 50 um layers or 
four 25 um layers. However, in the first case, a little less than 
twice as much time is spent in the freezing phase. Because of this 
the hailstone spends more time at 0°C and the total amount of heat 
transferred to the environment is somewhat larger when only two 50 um 
layers are accreted. This explains the slightly lower final hailstone 
averages temperature. 

However the main reason for the difference in the final 
surface temperature revealed in Table 5.8 can best be explained by 
considering Figure 5.9. When the first 50 um layer is accreted about 
2.4 J are conducted inside the hailstone in about .15 s. When the 
first 25 um layer is accreted about half of that amount of heat, 
roughly 1.3 J, is conducted inside during about one-quarter of the 
time (see Table 5.9). This means that, proportionately, the heat 
given to the hailstone in the last case is much closer to the surface, 
readily available to be transferred back to the atmosphere during 
the cooling phase (upper diagrams of Figure 5.9). When the 50 um 
layer is about halfway through the cooling phase, the surface is 
again raised to o°c in the other case, by the addition of another 
25 um layer (middle diagrams). Consequently there is more heat close 
to the surface, available for the next cooling period, and the surface 
temperature does not lower as much as for the 50 um layer (lower 


diagrams). In other words, the difference in final surface temperature 
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Figure 5.9 Hailstone temperature profiles during a 50 um growth. 
The left diagrams are for the accretion of one 50 um 
layer; the right diagrams two 25 um deposits. In 
both cases the initial conditions are a hailstone 
of 4.5 mm radius at -11°C in cloud condition a. 
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is due to the combined effects of raising the surface temperature to 
0°C more frequently with a larger number of small layers, and to the 
time lag required for the heat to diffuse back and forth inside the 
hailstone. 

The net effect of reducing (increasing) the thickness of 
the accreted deposit is to reduce (increase) the range of the surface 
temperature variations. Therefore returning to the hypothesis that 
the time cycle of the surface temperature can be viewed as the areal 
distribution of temperature over the hailstone, the effect should be 
reflected in the calculated standard deviation. In terms of such an 
interpretation the thickness of the deposit should correspond to the 
thickness of the droplets cr Cemethe | fF eacchet ion (droplets do spread 
upon accretion; Macklin and Payne, 1967) on the hailstone. Figure 5.10 
gives a plot of the standard deviation as a function of the thickness. 
As was shown in Figure 5.7, the standard deviations are slightly less 
under condition d, compared to condition a, indicating a weak depend- 
ence on the hailstone mean surface temperature. The two curves also 
indicate that the dependence on the accreted deposit thickness (or the 
size of the accreted droplets) is not very strong either, with a 
tendency to become even weaker for larger thicknesses. This is easily 
understood since the surface temperature cannot become lower than the 
air temperature nor higher than ote which puts a limit to the range 
of the temperature cycle. The reason for having the standard devia- 
tion go to zero is to indicate how this discrete model and the con- 
tinuous model could merge together, in the limit of an infinitesimal 


deposit. The limiting process could be visualized in the following 
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Deposit Thickness (10 m) 


Standard deviation of the surface temperature distribution 
as a function of deposit thickness (or droplet thickness 
after accretion) for cloud conditions a (upper curve) and 


b (lower curve). 
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manner. The infinitely small deposit (or accreted droplet) is warmed 
up to ac instantaneously, and cools back to the mean temperature in 
an infinitesimally short time. An additional calculation with a 
deposit thickness of 2.5 ym for condition d (plotted in Figure 2.10) 
seems to support this idea. 

Unfortunately, there is no simple way in which the above re- 
sults can be related to the size of the cloud droplets before accretion. 
Macklin measured experimentally the spreading of supercooled droplets 
freezing on an ice surface and found a very strong dependence on the 
temperature of the accreting surface (see Macklin and Payne, 1967), 
with the larger spreadings occurring at warmer temperatures. Other 
factors found to be of lesser importance were the speed of impaction, 
and the initial temperature of the droplets. 

Table 5.10 gives the droplet radius required to form a 
deposit of 25 um in thickness for a given mean surface temperature, 
as calculated by Macklin's method (Macklin and Payne, 1967). As can 
be seen droplets about three times larger are needed when the surface 
is about lace compared to -11°C. For the same mean temperature, 
when the other factors affecting droplet spreading are neglected, 
there is a linear relationship between the deposit thickness and the 
droplet radius before accretion. Therefore the conclusion to be 
tentatively drawn from Table 5.11 would be that a hailstone growing 
in conditions of large cloud droplets would have a "'broader'' surface 
temperature distribution than the same hailstone growing in a region 


where cloud droplets are smaller, other things being equal. 
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Table 5.10 


Droplet radius required to form a 25 um surface deposit, 
for different mean surface temperatures. Calculations 
were done using Macklin's method (Macklin and Payne, 


1967). 
aac) 02 0k 233) orc] -10.82 
r(um) 170 150 120 75 60 


Table 5.11 
Standard deviation o of the surface temperature 
distribution for different cloud droplets radii 


before accretion (calculated using Macklin's method). 
Calculations are for cloud conditions a and d. 


Cloud condition a 
o(°C) 1322 1.68 246 3.20 


r (um) 12 24 60 120 


Cloud condition d 


o (°C) 69 .89 et 1.47 T78 
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The effect of thickness was also investigated for the wet 
growth regime. No detectable influences appeared in the results as 
the calculated ice fractions were identical to the ones plotted in 
Figure 5.2, regardless of whether the deposit thickness was 5, 10, or 
20 um. This confirms the conclusion reached earlier in section 4.3 
that this discrete model becomes essentially similar to the continuous 


models when internal conduction becomes negligible. 


5.6 Freezing Times 


One of the fundamental assumptions of our model (after 
Macklin and Payne, 1967) is that the time of initial freezing is so 
short that it is negligibie compared to the time of subsequent freez- 
ing, and that no heat is exchanged outside the deposit during this 
first phase. As discussed in Chapter 1, the velocities for ice growth 
in freezing supercooled cloud droplets are not available; in order to 
obtain estimates of these velocities we use the linear crystalization 
velocities as a function of supercooling of Figure 1.1. Ratios of 
the estimated initial freezing time to the subsequent freezing time 
are given in Table 5.12. According to these values the initial 
assumption seems well justified. 

It might be more difficult to try to estimate the possible 
heat exchanges from the droplet during this phase. If we assume that 
the heat transfers due to forced ventilation calculated at 0° would 
constitute an upper limit (the dropiet should have an average tempera- 
ture lower than O°C during this phase), then the ratio of the heat 


transfers during the two phases would be tess than the ratios of the 
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Table 5.12 


Ratio of the estimated initial freezing time CreLoetiic 
subsequent freezing time tee one dia ferent. cloudeccn= 
ditions and mean surface equilibrium temperatures T . 
While this ratio is estimated to be an upper limit 

for the ratio of heat transferred to the environment, 
the square root of the ratio is estimated to be an 
upper limit for the heat conducted to the hailstone 

in both phases. 


ere 1 t/t, t/t, 
a -10.82 0040 062 
b =H ey 0025 050 
Cc =351| .0013 036 
d -2.08 .0010 031 
e -1.03 .0012 035 
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times. From Table 5.9 we find that the amount of heat conducted to 
the hailstone during the freezing phase is roughly proportional to the 
Square root of the time for which the surface is maintained at 0°C. 
Again it seems realistic that the extrapolation of this relationship 
to the initial freezing would yield an upper limit to the amount of 
heat conducted from the droplet to the underlying surface. Therefore 
the ratio of the heat conducted in both phases would be proportional 
to the square root of the time ratio, which is also given in Table 
5.12. Values are typically about 5 percent. It might be argued that 
this is not exactly negligible, thereby weakening the validity of our 
assumption. On the other hand, heat transfers outside the deposit 
are not simply disregarded. What our assumption does is to delay 
their realization till the subsequent freezing and assume that they 
occur at 0°C. But the effect of the delay is immaterial since the 
initial freezing time is neglected. Actually, what could happen, is 
that some of the heat would be conducted into the hailstone before 
the deposit reaches Gece thereby reducing the internal gradients and, 
perhaps, reducing the totai amount of heat conducted into the hail- 
stone. But certainly this effect would be much less than what is 


indicated in Table 5.12. 
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CHAPTER 6 


SUMMARY AND CONCLUSIONS 


Results show that even with the temperature cycle imposed by 
the model, the hailstone soon reaches a state of equilibrium in which 
the time average temperature is essentially equal to the equilibrium 
temperature as calculated by the continuous models. It is also found 
that the net effect of allowing internal conduction is negligible, 
yielding only a very slight average temperature gradient due to the 
slow rise of the surface average temperature of the growing hailstcne. 
However, within one accretion cycle there exists large heat fluxes 
of the same magnitude, but opposite in direction which have a de- 
terminating influence on the duration of the freezing phase and the 
overall time dependence of the surface temperature. Indeed, the 
surface temperature cycle adjusts in such a way as to try to balance 
the heat conducted to the hailstone during the freezing phase, with 
the heat transferred to the environment during the cooling phase. 

Under the hypothesis that individual droplets accreted by 
the hailstone would have a temperature cycle similar to the uniform 
deposit used in the model, an estimate of the temperature distribu- 


tion over the hailstone surface can be obtained. From such a 
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distribution, the mean surface temperature, the standard deviation and 
the fractional area of wet growth can be calculated. While the mean 
temperature is essentially equal to the equilibrium temperature of the 
continuous models, the standard deviation shows a weak dependence on 
the mean surface temperature and on the thickness of the droplets 
after accretion. From the fractional area of wet growth we see that 
the transition from the dry to the wet growth regimes occur fairly 
smoothly when the mean surface temperature rises towards Ogee and 
not in a discontinuous fashion, as implied by the continuous models. 

Although it is very difficult to analyze quantitatively 
the validity of the above hypothesis, qualitatively it can be argued 
that the overall effects of lateral heat transfers and lateral tempera- 
ture gradients would not change the freezing times or standard devia- 
tions to a very large extent. The stochastic nature of the accretion 
process has not been considered in this work. It is our feeling, 
however, that the possibility of overlapping and coalescence of accreted 
droplets over some area of the hailstone, combined with the fact that 
cloud droplets distribution is not monodisperse, would have a tendency 
to broaden the temperature distribution. 

As a final conclusion, this model] lends more support to 
the claim (Macklin and Payne, 1967; Pellet and Dennis, 1974) that the 
equilibrium temperature obtained from the continuous models is an 
adequate representation of the hailstone temperature. However this 
new model has the advantage of providing a simple estimate of the 


temperature distribution over the hailstone surface, which indicates 


that if the cloud droplets size distribution has little effect on the 
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mean surface temperature, it has some influence on the standard 


deviation. 
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20 
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11 


AaNANAN 


COMMON/BLOCK1/0 (101, 6) /ITENS/ETF, FFLUX,IR, IT, 
*R,TIM,ENMAX,TF,ENDIS(2) ,RINC 

REAL LWC 

DOUBLE PRECISION U 

READ(5,100) RINC,R,KMAX, IWRITE 

TR=100 

READ(5,200) TA,P, TO, LWC 

Tite et 

DR=(R+RINC) /IR 

DOMIOMK =A ET 

U (K, 1) =TO* (K- 1) *DR 

ie 2) St 

FUTIM=0. 

CWBAR=1. 0035+ .005*EXP (.693* (273-TA) /10.) 

EF(TASLT.253-)iCl2=. 235 

Te (DiGi 273 ed -6 07 

TRUTA. GE. 253. -AND-TA«LE.273.) C12=.207+.00140* 
* (TA-253.) 

PLUX=(1.44E-03% ((TA*P) **. 25) ¥(273-TA) +28.03% 
*C12* ( (TA** (.1)) J (P®*. 75) * (6. 107-VAP (TA) ))) 

RLF=79.7-CWBAR* (273.-TA) 


SoA eA CeR EEN GeeCrCl es 


DON le R= 4 KMAX 

DOM Jimena pated 

1 (N,5)=0. 

U(N,6) =0. 

R=R+RINC 

DR=R/IR 
VT=3705. *SQRT (TA¥R/P) 
TIM=3.668*RINC/ (VT*LWC) 
ELTIM=ELTIM¢TIM 

FIC FR=CWBAR*(273-TA) /79.7 
PNMAX=73.1*(1.-FICE) *RINC 
FFLOX=FLUX*.118/R¥**.25 
TF(IR. FQ. 1) GOTO150 
IF(U(II, 1) -GE.-- 1E-05) UCIT, 1) =-. 1E-05 
ETF=RLF*RINC*R/U (IT, 1) 
ETF=343. *ETF*ETP 

ro 


INITIATE FREEZING CYCLE 


CALL FREEZ (11,12,6150) 
IF(TF.GE.TIM) GOTO 59 
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G INITIATE COOLING CYCLE 


TC=TIM-TF 

IC=TC/TIM*100 

IF (IC.GT.50) Tc=50 

IF(IC.£Q.0)GoTO 59 

CAL ESCOOLM(LR iC; TC, TIN CID TAS De Retin) 


CALCULATE RESULTING ICE FRACTION 


@Q) @) OY & 


DFICE=1. 

GOTO 60 
59 DF ICE=FICE+ (1.-FICE) *ENDIS (12) /ENMAX 
60 CONTINUE 

IF (MOD(K, IWRITE).GT.0)GOTO 7 


OG e UO Baek eS ULL 


EG Raina 


WRITE (6, 300) TIM, ELTIM,DFICE, ETF 
RR=0. 

DOMCMN=o ite 2 

RR=(N-1) *DR 

TEMP=U(N,1I2) /RR 

IF (IC. FO. 0) U(N, 4) =0. 

TEMP1= (N,4) /RR 

ATEMP=(U (N, 5) ¥TF+0 (N, 6) ¥TC) /TIM/RR 
TEMPF=0 (N,3) /BR 

WRITE (6,400) RR, TEMPF, TEMP1, TEMP, ATEMP 
CONTINUE 

IF (K. NE. KMAX) WRITE(6,600) 

CONTINUE 


(oe) 


REINITIALIZATION OF THE GRID AFTER ACCRETION 
OF ANOTHER LAYER 


QEOI@Oi@i@u— 


DO 9 N=2,IR 
9 U(N,11) =(U(N+1, 12) -0 (N-1,12) ) *RINC*(N-1)/ 
*(2.*IR¥*DR) +0 (N, 12) 
U(IT,11) =U (IL, 12) *(R+RINC) /R 
TF(i1.80.1) GOTO 1 
NOed2eNa lei. 
12 U0 (N,1) =U (N,11) 
GOTO 1 
459 DFICE=FICE+(1.-FICE) *FFLUX*TIM/ENMAX 
WRITE (6, 300) TIM, ELTIM, DFICE 
1 CONTINUE 
STOP 
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100 
200 
300 


400 
600 


PORMAT (2F10.8,415) 

PORMAT(3F10.5,P710.9) 

FORMAT (1H ,"FALL TIME',F7.2,5X,*ELAPSED TIME? 
*,F7.2,5X%,"DEPOSIT ICE PRACTION',F5.2,F12.-4) 
PORMAT (1H ,1F10.5,4F11. 3) 

FORMAT (1H1) 

END 
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THIS FUNCTION CALCULATES THE SATURATION VAPOOR 
PRESSURE WITH RES PECT® TOW WATER, ACCORDING TO 
WEG LEV Seis THODMO AN, a3 ,ee ee UG. 


FUNCTION VAP(T) 

TH=1.-373/T 

X1=TH 

X2=X1*TH 

X3=X2*TH 

A=13.3185*X1-1.976¥*X2-. 6445*X3-.1299*X3*TH 
VAP=1013.25¥*EXP (A) 

RETURN 

END 
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THIS SUBROUTINE SOLVES THE HEAT EQUATION WITH 
THE SURFACE TEMPERATURE MAINTAINED AT OC UNTIL 
THE SURFACE DEPOSIT IS ENTIRELY FROZEN. 

THE CRANK-NICHOLSON METHOD IS USED. 


SUBROUTINE FREEZ (11,12,*) 
COMMON/BLOCK 1/U (101, 6) /ITEMS/ETF,FFLUX,IR, 


SIT Rep LLM, ENMAX, DFO ENDIS (2)), RING 


DIMENSION CSTAR (500) ,DSTAR (500) 
DOUBLE PRECISION U,CSTAR, DSTAR, ALM,GAM,DR2, 


BUCO GT 


DR=R/TR 

TH=.5*IT 

DT=ETF/IT 
IP(ETF.GE.TIM) DT=TIM/IT 
TPCETP.GE.TIM/2.) LH=1 
ALM=.0109*DT/ (DR*¥DR) 
BETA=ALM/(ALM+2.) 
ALM=ALM*(1.-1./(2.+ALM) ) 
ENDIS (1) =0. 

ENDIS (2) =0. 


CALCULATE THE CINITIAL HEAT CONTENT 


DR2=DR*DR 

HCO=0. 

90.122 SN=2/ 20k 
HCO=HCO+ (N-1) *0 (1,1) 
HCO=(HCO+IR¥*U (IR+1,1) /2.) *DR2 
HCO=HCO-RINC#*U(IR+1,1) *R 

IF (HCO/(4.189*R*R¥*R) .GE.-.001) GOTO 150 
HCO=.45*HCO/ (R*B) 
GAM=-ALM*¥2.-1. 

MLIM=2*IT 

CSTAR (1) =ALM/GAM 


U(TR+1,1) =0. 
FK=RINC/DRt. 2 
IK=FK 

IF (IK.EQ.0)GOTO 45 
HO rsoms— lk 
M(TR+ iJ 91} -0- 
CONTINUE 


Py ib date 
sea” 


Ad 


“er onte tro 


ft ee 


QiGu@ 


AA Gea Q @) 


Wey QE 


we) 


MR AQ @ 


104 


START TIME ITERATIONS 


DO1 M=1,MLIM 
TF=M*DT 
I2=MOD(M,2)+1 
mia? 

IF (1I1.£0.0)11=2 


CALCULATE THE TRANSFORMED COEFFICIENTS BY THE 
THOMAS ALGORITHM 


DST AR (p= (U (2,01) (1-25 *B ETA) 4B ETA+0 (3, 04))) 
*/GAM 

IT=IR-1 

DOM a2 tL 

CSTAR (I) =ALM/ (GAM-ALM*CSTAR (I-1)) 

DSUAR Gia (=U 17 1) li. - 2. BETA) 
SEB PAS (UCL Ol) tt (Lt2, 0.7) 11) 
% APN DS CAR P13) 7 (GAN-ALM*CS DNR T= 130) 


CALCULATE VALUES AT GRID POINTS 


U(LPR,12) =DSTAR IT) 
U¢IR+1,12) =0. 
U(URyoON=U (TR, 12) t0 (18,5) 
U(IR+1,5) =U (IR+1, 12) +U (IR+1,5) 

DO3 N=2,I1 | 
U(IR+1-N,12) =DSTAR(IR-N) -CSTAR(IR-N) * 

TU CLEP Z—N ys) 

Uc IR+1-N, 5) =U(IR+1-N,5) +U( IR+1-N, 12) 
IF(M.LE.IH)GOTO 1 


CALCULATE RESULTING HEAT CONTENT 


HCT=0. 

DO 33 N=2,IR 
HCT=HCT+ (N-1) *U (N, 12) 
HCT=(HCT+IR*U (I R+1,12) /2.) *DR2 
HCT=.45*HCT/ (R*¥R) 
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CALCULATE DISSIPATED ENERGY 


Qa@ie 


ENDIS (12) =FFLUX*¥M*DT+HCTxHCO 
c 
c CHECKERED EPOS LETS PENTIR EGY sPROZEN 
C 


IP (ENDIS(1I2).GE.ENMAX)GOTO 11 
IF(TF.GE.TIM) GOTO 6 


1 CONTINUE 
11 CONTINUE 
TP ( (ENDIS (12) -ENMAX) .GE. (ENMAX-ENDIS (T1))) 
*T2=71 
6 CONTINUE 


TF=TF-.5*DT 
OF LUX=PFLUX*TF 
L=IR+1 
DOS el=—1i0L 
Pree. 2) U (1, 1) =U (1, 2) 
Ot) =U (1,5) /i 
qh Ut oy Suis) 
WRITE (6, 100) TF, DT, ENMAX ,ENDIS(1I2) ,OFLUX,HCO, 
=HOd 
100 BORMATNO HOP SE 12.5, 2020.0) 
RETURN 
150 TR=1 
RETURN1 
END 


or * 2 
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THIS SUBROUTINE SOLVES THE HEAT EQUATION WITH THE 
SECOND BOUNDARY CONDITION. 
THE CRANK-NICHOLSON METHOD IS USED. 


LENE ES Is Selene. Aue itd. Wey Wee oe ae 
{2) 

COMMON/BLOCK 1/0 (101, 6) 

DIMENSTON CSTAR (100) ,DSTAR(100) ,TT(50) 
DOUBLE PRECISION U,CSTAR,DSTAR,ALM,GAM,HCF, 
*DR2,D1 

DT=TC/IT 

THALF=1T/2 

DR=R/IR | 

ALM=.0109¥*DT/ (DR¥*DR) 

BETA=ALM/(ALM+2.) 

ALM=ALM* (1.-1./ (ALM+2.)) 

GAM=-ALM*2.-1. 

VFLUX=0. 

MLIM=IT 

CSTAR (1) =ALM/GAM 

Q1=3. 39F-02* ( (T*¥P/R) **. 25) 
02=660.*C12*(T#*. 1) /( (DEX. 75) ¥ (R25) ) 
ESV=VAP(T) 


START ‘TIME ITERATIONS 


DO 1 M=1,MLIM 
T2=MOD (4,2) +1 
T1=12-1 

Tec ee OO} ct =2 


GALCULATE THE TRANSFORMED COEFFICIENTS BY THE 
THOMAS ALGORITHM 


DSTAR (1) =- (U (2,11) * (1.-2. *BETA) #BETA*U (3,11) 
*) /GAM 

TI=IR-2 

n@) 2 aR ime 

CSTAR (I) =ALM/ (GAM-ALM*CSTAR (1-1) ) 
DI=-(U(I+1,11)* (1.-2-*BETA) #+BETA* (0 (1,11) + 
*0 (I+2,11))) 

DSTAR (1) =(DI-ALM*DSTAR (I-1)) /(GAN-ALM* 
*CSTAR(I-1)) 
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CALCULATE VALUES AT THE TWO INNER SURFACE 
POINTS BY THE METHOD OF FALSE POSITION 


DI=-(U(IR,11) *(1.-2.*BETA) +BETA* (U(IR-1,11) 
*+U (IR+1,11))) 

A= ((DI-ALM*DSTAR(IT)) /(GAM-ALM*CSTAR (IT)) 
*—-DSTAR (11) /(4.*+CSTAR (IL) ))/8 
B=(3.-2.*DR/R) /(4.+CSTAR (II) ) +ALM/(GAM-ALM*® 
*CSTAR(IT)) 

FXO=DR¥* (Q1*(273.-T) +02* (6. 107-ESV)) *2./(4.+ 
*CSTAR(II))-A 

XN=T-273. 

Domai — 120 
ESVI=6.109177956+XN* (5.03469 897E-014+XN* 
*(1.886013408F-02+4XN* (4. 176223716EF-04 +XN* 
*(5.82472028F-06+XN* (4. 838803174E-08+ 
*XN*1.838826904F-10))))j) 

PXN=B*XN+DR* (O1* (XN4+273.-T) +02* (ESVI-ESV)) 
**2./(4.4+CSTAR(TI))-A 

XN1=-XN*FXO/( FXN-FXO) 

AV=XN1-XN 

Tr (aA¥.Le..004j)GOTO 4 

XN=XN1 

CONTINUE 

U(IR+1,12) =XN1*R 

ESVI=6.109177956+XN 1% (5.03469897E-01+4XN1* 
*(1.8860134 08F-024+XN1* (4.1762237 16 E-O4+KN1* 
*(5.82472028E-06+XN 1* (4. 838803174E-08+ 
*XN1*1.838826904E-10))})) 

OV=O1* (XN14273. -T) +Q2*(ESVI-ESY) 

0 (IR, 12) = (DI-ALM* (0 (IR+1,12) +DSTAR{II)))/ 
* (GAM-ALM*CSTAR(IT) ) 

VFLUX=VFLUX+. 0O05S*QV*DT 

OU (IR+1,6) =U(TR+1,6) +U (IR+1,12) 

U(IR, 6) =U(IR,6) +U(IR,12) 

II=IR-1 


CALCULATE VALUES AT INTERIOR POINTS 


pO 5 N=2,II 

0 (IR+1-N,12) =DSTAR(IR-N) -CSTAR(IR-N) * 
*U (IR+2-N,12) 

U(IR+1-N, 6) =U (IR+1-N, 6) +U (IR+1-N, 12) 

PT (My) =XN1 

IF(M.NE.IHALF) GOTO 1 
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e STORE VALUES AT THE TENTH TIME STEP 
c 
IN=IR+1 
DO 66 N=1,IN 
66 U(N,4) =U(N,1I2) 
1 CONTINUE 
re 


DO 99 N=1,IN 
g U(N,6) =0 (N,6) /IT 


=) 
if 
c CALCULATE FINAL HEAT CONTENT 
e 


DR2=DR*DR 
HCF=0. 
DO simN= 2 ER 

Ri HCF=HCF+ (N-1) #0 (N,12) 
HC F=(HCF+ (IR) *U (IR+1,12)/2.) *DR2 
HCF=HCF*.45/ (R¥R) 
WRITE (6, 100) TC, VFLUX,HCF,AV 
WRITE (7,200) (TT (1) ,l=1,17) 
Tie=0- 
TB2=0. 
DO 88 I=1,1IT 
TB=TB+TT (I) 

88 TBR2=TB2tTT (I) *TT (I) 
TB=TC/TIM¥TB/IT 
TB2=SORT (TC/TIM*TB2/IT-TB*TB) 
WRITE (7,260) TB, TB2 

200 FORMAT(5X,15F8.2) 

400 FORMAT (1H0, 2P10.5,020.5,F30. 10) 
RETURN 
END 


APPENDIX 2 


109 


110 


LI STOGRAM 


The Listogram gives solutions of equation (1.2) - actually, 
this equation divided by we - in graphical form. For a given air 
temperature and liquid water content, the hailstone temperature and 
ice fraction are found or interpolated from the labeled curves on the 
graph. The air temperature, pressure and height are related accord- 
ing to a simple cloud model given in List et al. (1965). Therefore 
one must be careful to select the appropriate temperature and pressure 
in order to calculate the Reynold number Re, (equation (2.13)). How- 


ever, List et al. (1967) claim that 


“lL Re = = 55.84 0-75 


with an accuracy of better than + 13 percent. The product Wed allows 

to get an easy estimate of the hailstone temperature from the diagram. 
The heavy continuous lines on the Listogram define the 

regions where one heat transfer is dominant. The ratios are defined 


as follows: 
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The dashed lines define regions where those transfers are 


absolutely dominant, that is, greater than .5. 
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